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I. INTRODUCTION 

This report summarizes the work done on the extension of Contract 
NAS8-28019 during the period of January - December, 1973. The presen- 
tation here continues that of the Interim Report [1]* on the same contract 
published May, 1973, Much of the background for this report is contained 
in the Interim Report. 

Briefly, the Interim Report presented feasibility results for 
smoothing surface irregularities of a homogeneous mirror by using a con- 
trolled pattern of heating units on the rear of the mirror. It was shown 
that significant deflections of the mirror surface could be caused. How- 
ever, no actual error-suppression simulations were made. 

This report presents the generalization that was made to model a 
layered structure of a kind that represents a light-weighted mirror. This 
theory is contained in Chapter II, In addition, the strategy for error 
suppression is derived in Chapter III. The results of a variety of error- 
suppression studies are presented in Chapter IV. The computer programs 
for all parts of this study are contained in Appendix A and Appendix B. 


* Numbers in brackets indicate entries in REFERENCES. 




II. THERMAL RESPONSE ANALYSIS 
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The thermal response model of the primary mirror is similar to the 
model that was developed and presented in the Interim Report [1] of this 
project. The model was generalized to allow the primary mirror to be a 
three-layer structure as sketched in Figure II-l. The thermal and me- 
chanical properties of each layer are orthotropic. The three orthotropic 
layers of the primary mirror represent the front, core and back. The 
orthotropy of each layer (particularly the core and back) is included to 
account for the directionality of the equivalent mechanical and thermal 
properties induced by the light-weighting of the mirror. 



Figure II-l Layered Model of Primary Mirror 
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Thermal Properties of Each Layer 

The thermal properties of the i-th layer (i-1, 2, 3) are the 
following: 

C - Specific thermal capacity which is equal to the spe- 

cific heat times the mass density, 
k . k_, k, - Thermal conductivities in the r, 0, z- directions. 
a r’ a e* a z " Coe ff lc i ents of linear thermal expansion in the 
r, e, z- directions. 

These thermal properties are constant within each layer. 


Mechanical Properties of Each Layer 

The mechanical properties of each orthotropic layer are given as 
the stress-strain relations given in [2], 


e r = ( 1 /E r )o r * (v zr /E z )o z 

E e = ' (v ra /E r )0 r + ( 1 /E e )o e ' (v ze /E z )o z 
e z = ' (v rz /E r )a r ' (v ez /E e )a a + ( 1 /E z )o z 


(II-D 


where 


re 


rz 


ez 


( 1 /^re) T re 
- <’' G rz>Vz 
■ ^ezKz 


a r > ° 0 , 5 Z 5 normal stresses 

T re* T rz » x q z = sh ear stresses 

e = normal strains 
r b z 

Y ro ’ Y rz 5 Y 0z 


( I 1-2 ) 


= shear strains 
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and the elastic properties E . E_, E , must satisfy the 

r y z or uz 

relations 


E v 
r er 


= E a u ra ’ E g v ze 


= E , E v 
z ez # z rz 


= E r v zr 


(n-3) 


The Discrete Model of the Primary Mirror 

The temperature and the displacement components are expressed in 
Fourier series forms. The Fourier coefficients of these four fields 
vary over the meridional (r,z) section. The (r,z) variation is dis- 
cretized by the finite element method. This method of discretizing the 
temperature T(r,e,z,t) and the displacements u (r,e,z), u A (r,e,z) and 
u 2 (r s e,z) is exactly the same as previously used in our Interim Report 
[1], The expressions for the temperature and the displacement fields 
are 

N 

T(r,e,z,t) * ^T n (r,z,t)e 12mT0 (II-4) 

n=-N 
N 

u r (r t e,z,t) * ^U n (r,z,t)e l2n7T0 (II-5) 

n=-N 

N 

u Q (r,0,z,t) = (r,z 5 t)e l2mT0 (II-6) 

n=-N 

N 

u z (r,0,z,t) = ^]w n (r,z,t)e l2n7re (II-7) 

n=-N 


where the component fields T p , U n , V n> W n are such that T, u r , u Q and 
u z are real valued functions. 

Each of the component fields is discretized over the meridional 
section by assuming linear variation over arbitrary triangular areas. 
The section is divided into a collection of triangles as sketched below 
in Figure II-2. 
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Figure 1 1-2 Meridional Section of the Mirror 
Consider any one of the triangular areas as sketched below. 



Figure I 1-3 Triangular Area 


The component temperature T n (r,z,t) is represented as 

T n (r,z,t) = [*,(r,z),<f. 2 (r,z),<t> 3 (r,z)]T fl (t) { 1 1-8) 

where T^t) is a matrix of the nodal values of the component temperature 
yt> = [T„i (t) , T n2 (tl,T n3 (t)] t (11-9) 
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and [c] t indicates the transpose of [c]. The interpolating functions 
^ (r ,z) are given by 

(r,z) - a i + b^r + c^z (11-10) 

where 


a 

b 


c 

A 


1 

1 

1 

t 


( r 2 z 3 ' r 3 z 2^ A t 


< z 2 - Z 3 )/A t 
(r 3 - r 2 )/A t 

[( r 2 - r i)( z 3 - z p 


(i" 3 - r,)(z 2 - z^)V 2 


(11-11) 


( 11 - 12 ) 


for i = 2 and 3 permute the indices. 

Likewise the component displacements are represented as 
U n (r,z,t)1 

4 V n (r,z,t)>= [<J»-| (r,z)I, <j» 2 (r,z)I, <|> 3 (r,z)|] U^t) (11-13) 

k W n (r,z,t)j 

where I is a 3x3 identity matrix and U (t) is a matrix of the nodal 

-n 

values of the component displacements given by 

V tJ = CU nl (t) , V nl (t), W nl (t). U n2 (t) W n3 (t)] t (11-14) 

Utilizing the discretizations of equations ( I 1-8) and (11-13) over all 
the triangular areas that comprise the meridional section (r,z - section) 
of the mirror, the component temperature and displacements can be char- 
acterized by the values they have at the nodes of the section. 

The temperature and displacements are governed by sets of 
partial differential equations and appropriate boundary conditions. 
Alternatively, the temperature and displacements are formulated by vari- 
ational principles that are generalizations of those given in the In- 
terim Report [1], The generalizations involve accounting for the 



orthotropic nature of the layers of the mirror. The development dif- 
fers so little from that of the Interim Report that it will not be re- 
peated here. 

The finite element model of the component temperature results in 
a set of matrix equations of the form 

^ In + B t, = % Ul-15) 

for each Fourier component. Similarly, the component displacements are 
governed by 

K U + P T = 0 (11-16) 

-n -n -n -n ' ' 

Thus equation (11-15) is integrated for each set of nodal components 
T^(t) which are substituted into equations (11-16). Equations (11-16) 
are solved by matrix inversion for U^(t). The nodal components are sub- 
stituted into equations ( I 1-8) and (11-13) to obtain T n (r,z,t), U n (r,z,t), 
V n (r,z,t), and W n (r,z,t) which are then substituted into equations ( I 1-4) , 

( 1 1-5) , (1 1-6) and ( 1 1 -7) to yield T(r,e ,z,t) , u (r,e,z,t), u fl (r,6,z,t), 

r y 

and u z (r,e,z,t). 

The component heat input comes from the heating that is applied 
to the rear of the mirror. To characterize the steady-state thermo- 
elastic response of the mirror, the normal displacements at a set of 
points of the front surface of the mirror due to placing a unit heater 
over a patch of the back were calculated. 

The points on the front surface (the sample set) are at the inter- 
sections of the circles and radial lines shown in Figure 1 1-4. The 
heater patches are between the concentric circles and are centered 
about the radial lines as shown in Figure 1 1-5. The indices of a heater 
is the index of the inner circle and of the radial line. Let I be the 
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• - Support Points 


Figure II-4 Sample Point Locations 
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index of the concentric circles (I = 1 is the inner boundary and I = IM 
is the outer boundary). The radial lines are indexed by K (K=l is for 
6=0, K=KM for e~ ( KM- 1 )/KM*360°) . The sample points are numbered from 
1 to IM x KM. A particular sample point is identified by the index ISP 
given by 

ISP = (K-l)IM + I I = 1,2,..., IM (11-17) 

K * 1,2,..., KM 

Likewise the heater patches are numbered from 1 to (IM-l)xKM. The index 
IHP of a heater is given by 

IHP = (I-l)KM + K I = 1 ,2 IM-1 (11-18) 

K = 1,2,... ,KM 


The Thermoelastic Influence Matrix 

Define the steady-state normal displacement of the front surface 

at i due to a unit heat patch at j as a. .. The displacements w at the 

' J 

sample points due to heat inputs ^ at the heater patches are given by 
the matrix equation 

w = Aa (II-19) 

For the calculations carried out in this study, IM = 15 and KM = 12 so 
that there were 180 sample points and 168 heater patches. This in- 
fluence matrix ft represents the statical thermoelastic response of the 


mi rror. 



III. ANALYSIS OF ERROR SUPPRESSION 
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Introduction 

The influence matrix which is computed from the structural response 
program may be used to associate a controlled deflection of the sample 
point set to a set of heat inputs on the back face of the mirror. The 
purpose of the control program is to compute the amount of heat input 
to a designated pattern of heat patches that reduces an arbitrary sur- 
face error to its minimum in an r.m.s. sense. The control program also 
computes the surface error (r.m.s.) before and after compensation. The 
r.m.s. error is interpreted as the root mean squared difference of 
sample point location of the mirror surface and the best-fit-sphere of 
the surface passing through the rigid support point. 

Removal of Best Fit Sphere 

In heating the mirror for control purposes, a significant amount ‘ 
of heat goes into the free thermal expansion of the mirror. This de- 
flection does not contribute to local error smoothing and may be com- 
pensated by refocusing of the mirror. Moreover, if retained as part 
of the control, it significantly reduces the sensitivity of the con- 
trol function. For these reasons the free thermal expansion terms are 
measured and their effect is removed from the influence matris. 

Consider the mirror geometry shown in Figure III-l, which shows a 
spherical surface of radius R t that is the best fit sphere through a 
set of displacements w_ measured from the perfect sphere of radius R 
and also passes through the support points of the mirror. 
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Figure III-l Geometry of Best-Fit-Sphere 
The difference in the two reference spheres is 

y = < R i- R > 

The best fit sphere determined by the radius R-j- is found by minimizing 

the functional T , 

IM oM 

j= Z&ij-v 2 

i=l j=l 

1. L 

where w. . is the displacement of the ij— surface nodal point due to 
i J 

heat input on the rear and y^ is the difference at e... 
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The minimization yields 

IM IM , IM 

c-* r-« COSO- / r-i COSO. 2 

A = R T - R = 2_ 2_ "ijfenr - V - U' 

1-1 1-1 0 ' i-i 0 

which represents the adjustment of the reference sphere. 

The influence matrix as computed in the response program represents 
the correspondence of sample point deflection to heat input. Since the 
mirror may be focused to remove spherical errors, that part of the de- 
flection which may be considered a readjustment of the reference sphere 
may be ignored. Therefore each column of the influence matrix (corre- 
sponding to one heat patch) may have its best-fit-sphere of influence 
removed. This is performed on each column of the influence matrix in 
the control program by application of (III-l). The modified influence 
matrix represents the influence of heater patches on sample point de- 
flections about the resulting best-fit-sphere. 

When the sample point uncompensated error vector is computed for 
the particular error under consideration, its best-fit-sphere is simi- 
larly removed. This corresponds to re-focusing the mirror prior to 
error suppression. The process is carried out in the control program. 

In addition, the control program reduces the modified influence 
matrix to the appropriate array for the sample set and heat patch set 
selected. 

Error Computation and Suppression 

The deflection of the front of the mirror about its best-fit- 
sphere is composed of a disturbance term w^, whose effect is to be 
minimized and a controlled deflection term imposed to reduce the 
effect of Both w's represent a vector of sample point displacements 
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In the control program the modified and reduced influence matrix & 

is used to associate the heat patch input vector to the controlled 

deflection w . 

-c 

«c = AS- 

The total deflection vector is 
w = 

The purpose of the heater pattern is to produce a vector £ that 
minimizes the r.m.s. error 

J = /(wV)/N 

where N is the number of sample points. It is understood that the modi- 
fied influence matrix and modified disturbance vector are used so that 
the resulting error is measured about the best-fit-sphere. 

The minimization of the r.m.s. error is accomplished by 

a = - (aV 

This distribution of heats on the rear of the mirror yields a com- 
pensated r.m.s. error of 

j c - 4^0 - ataV a^v" 

or a reduction in J 2 of 

J2 ■ J c 2 = "d^^ar 1 aV N 

The resulting surface error figure for a particular disturbance 
was normalized by the uncompensated error for the same disturbance. 

The resulting r.m.s. error figure is 


15 


J = . 1 


*4 ^ 


Heat Input and Deflection Measurement Points 

The node distribution for both the front and the rear of the mirror 
is shown in Figure I 1-4 and II-5. The surface has IM concentric nodal 
rings and KM nodal rays. The nodes occur at intersections of rings and 
rays. The nodes are selected by the program after mirror dimensions, 

KM, and IM are given. 

The nodes for measurement of surface deflection are chosen in the 
control program. Any subset of the full set of nodes may be chosen. 

The heat is applied to the rear by patches. The patches extend 
the entire distance between rings and make an angle PAN centered about 
a ray. Patch location and patch angle may be inputed in the control 
program. 

To examine the feasibility of error suppression, the maximum set 
of sample points was used in the runs. Heater distributions comprising 
from 6 to 48 heat patches were used. The arrays of the heat patches 
are shown in Figures III-2 through I I 1-6 . 





























IV. PRINCIPLE FEASIBILITY RESULTS 


Introduction 

The influence matrix corresponding to the mirror described in 
Table IV-1 was computed using the Response program outlined in Chap- 
ter II and Appendix A. The Control program described in Chapter III 
and Appendix B was used to compute the heat input pattern, unsuppressed 
surface error, and suppressed (controlled) surface error for a variety 
of surface errors which are described in the following section. For 
each error, patterns of 6, 12, 24, 36, and 48 symmetrically located 
heaters (see Figures I I 1-2 through I I I -6 ) were used. The tabulation of 
results is given in Table IV-2 and discussed in a later section. In 
each case the error corresponds to the r.m.s. deviation of 180 sym- 
metrically located sample points, (see Figure I I -4) , about the best- 
fit-sphere to the points that passes through the supports. 

Uncompensated Surface Errors 

Five typical surface errors, which were suggested by M.S.F.C., 
were used for feasibility studies. These are shown in Figures IV-1 
through IV-5. The descriptive names of the disturbances are: 

(A) High in the middle 

(B) Warp 

(C) High ring 

(D) High spot near support 

(E) High spot between supports 

Error-Suppression Results 

The results obtained by using the control program to suppress the 

five errors are shown in Table IV-2. A set of 180 unweighted sample 

21 



Table IV-1 


Dimensions and Properties of NASA Test Mirror 


Dimensions : 

Inner Radius 
Outer Radius 
Thickness 
F-number 


0.0 

10.0 in. 
0.667 in. 
2.0 


* Properties : 

Specific Heat 

Conduct! vi ty 

Emissivity 

Young's Modulus 

Poisson's Ratio 

Coefficient of 

Thermal Expansion 


0.015 Btu/in 3 - °F 
0.05 Btu/br - in - °F 
0.04 

1.06 x 10 8 lb/ in 2 
0.17 

0.311 x IO'Vf 


* This mirror is homogeneous and uniform so the properties do not 
depend upon direction and are constant throughout the mirror. 
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Support Points 



Figure IV-1 Disturbance A - High in the Middle 
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•- Support Points 


w d (x) = W Q (r s - x)(x + r $ ) 


Figure IV-2 Disturbance B - Warp 
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Figure IV-3 Disturbance C - High Ring 



w d (x,y) = 


• - Support Points 


- [(X - r s /2)2 + y2] 1/2 


/ 


I 0 ^ j. r s /4 

w [1 + cos(^™)], r-j < r $ /4 

L S 


Figure IV-4 Disturbance D - High Spot Near Support 



Number 


Fractional Error with Disturbance 


Heat Patches 

A 

B 

C 



D 

E 

6 

.995 

.734 

.587 

.891 

.876 

12 

.992 

.551 

.576 

.842 

.846 

24 

.991 

.481 

.362 

.801 

.811 

36 

.991 

.449 

.278 

.729 

.753 

48 

.991 

.434 

.268 

.562 

.580 
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points was used. The five heat patterns described earlier were used 
for each error. The error of the compensated mirror is, in each case, 
normalized by the uncompensated error. Each error is with respect to 
the best-fit-sphere as described earlier. 

A substantial improvement occurs for all errors but the first. The 
effect of increasing the number of heaters differs from one error to 
the next. Generally the more localized errors responded well for the 
larger number of heaters and not particularly well for the small heater 
patterns . 

The less localized errors of warping and high ring were relatively 
correctable with a small number of heaters. The flattened mirror (A) 
did not show any appreciable degree of correctability with any of the 
heat patch distributions. The error reductions as a function of the 
number of heat patches for the disturbances are presented in Figure IV-6. 
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APPENDIX A 

DESCRIPTION OF THERMOELASTIC RESPONSE PROGRAM 

Purpose 

This program computes the thermoelastic deflection of a set of nodes 
specified in an axi symmetric body. Heat is inputed on one surface of the 
body and radiates at the other. 

Options 

This program provides the following options: 

(1) The temperature and deflection of the node set may be computed 
at any time. 

(2) The steady-state temperature and deflection may be computed. 

(3) The influence matrix of surface deflection for a designated 
heat input pattern may be computed. 

Input Parameter Definition 

Definition 

The total number data cases to be run. 

Number of node locations in the radial 
direction. (Cannot exceed 15.) 

Number of node locations in the longi- 
tudinal direction. (Cannot exceed 15.) 
Radiation coefficient of the front surface. 
Temperature of radiating reference medium. 
Integration time. 

Number of integration steps. 

Number of integration steps between 
printings. 


Parameter 

NUMPAT 

IM 

JM 

CF 

TO 

DELT 

NTS 

IPRINT 
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Parameter 

TI 

NM 

IS 

KM 

52 

53 

INFLU 


IP 

KP 


Defi ni tion 

Initial mirror temperature. 

Highest harmonic (angular) component. 
Location of radial node of support ring. 
Number of angular divisions. Must not 
exceed 12. 

Angular position of second support (first 
one is at 0). 

Angular position of third support. 

Switch to choose simulation of response 
or influence matrix computation. 

If INFLU = 0 Program calculates the re- 
sponse to heating one patch. 

If INFLU > 0 Program calculates the 
thermoelastic influence coefficient 
matrix. This matrix is written in a 
data file through 1-0 unit 4. 

Radial node of single patch input. 

Angular position of single patch input. 
Patch angle for single patch. 


PA 
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Parameter 

Definition 

PH 

Heat input rate for single patch. (Negative 
is input; positive is output.) 

PAN(I) 

Patch angle for i-th radial node in full pat- 
tern generation. 

DO 

Outside diameter of mirror. 

DI 

Inside diameter of mirror. 

H 

Thickness of mirror. 

FNO 

F-number of the mirror= Focal Length/Diameter 
of Mirror. 

E(1,M) 

of layer M (M=Bottom, Core, Front) 

E(2,M) 

E 0 of layer M. (See equations II-l) 

E(3,M) 

E z of layer M. 

G(1,M) 

of layer M. 
ro J 

G(2,M) 

G _ of layer M. 
rz 

G(3,M) 

G of layer M. 
0 z 

V { 1 ,M) 

v 0r of layer M. 

V(2,M) 

v zr of layer M. 

V(3,M) 

v of layer M. 

ALF(l.M) 

of layer M. 

ALF(2,M) 

a of layer M. 

ALF(3,M) 

a z of layer M. 

CK(] ,M) 

k^ of layer M. 

CK(2 ,M) 

k of layer M. 

0 

CK(3,M) 

k z of layer M. 

CP(M) 

Specific thermal capacity of layer M. 
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Input Data Card Listing 


Card No. 

Parameter 

Data Field 

Format 

1 

NUMPAT 

1-5 

15 

2 

IM 

1-5 

15 


JM 

6-10 

15 


CF 

11-20 

F10.5 


TO 

21-30 

F10.5 


DELT 

31-40 

F10.5 


NTS 

41-45 

15 


IPRINT 

46-50 

. 15 


TI 

51-60 

F10. 5 

3 

m 

1-5 

15 


IS 

6-10 

15 


KM 

11-15 

15 


S2 

16-20 

15 


S3 

21-25 

15 


INFLU 

26-30 

15 

4 

IP 

1-5 

15 


KP 

6-10 

15 


PA 

11-20 

F10 


PH 

21-30 

F10 

4 

PAN ( 1 ) . . . PAN ( 7 ) 

1-70 

7F10.5 

5 

PAN ( I M- 7 ) PAN ( I M- 1 ) 

1-70 

7F10.5 

6 

DO 

1-10 

F10. 5 


DI 

11-20 

F10.5 


H 

21-30 

F10. 5 


FNO 

31-40 

F10.5 



The following cards contain the material properties for each layer. 
There should be three sets of these two cards. 


E(1,M) 

1-8 

F8.5 

E(2,M) 

9-16 

F8.5 

E(3,M) 

17-24 

F8.5 

G(1,H) 

25-32 

F8.5 

G(2,M) 

33-40 

F8.5 

G(3,M) 

41-48 

F8.5 

V(1 ,M) 

49-56 

F8.5 

V(2,M) 

57-64 

F8.5 

V(3,M) 

65-72 

F8.5 

ALF(1,M) 

1-8 

F8.5 

ALF(2,M) 

9-16 

F8.5 

ALF(3,H) 

17-24 

F8.5 

CK(1 ,M) 

25-32 

F8.5 

CK(2,M) 

33-40 

F8.5 

CK(3,M) 

41-48 

F8.5 

CP(M) 

49-56 

F8.5 


Card 2 through the end constitute one data set. There should be 


NUMPAT data sets. 



Output of Program 
A. Option 


INFLU = 0 
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1. Repeated input data. 

2. Table of temperatures and deflection of the node (I=IM, O-JM) 
for each node for a check of convergence. 

3. Table of displacements and temperatures of the nodes on the 
front surface before the support constraints are added. 

4. Table of displacements of the nodes on the front surface after 
supports are added. 

B. Option INFLU = 1 

1. Repeated input data. 

2. Patch angles of each patch ring. 

3. The displacements of the nodal points on the front surface of 
the mirror corresponding to each patch location are read into 
a f i 1 e . 

Program Listing 

A listing of the program and its subroutines appears on the 


following pages. 
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C PROGRAM TO DETERMINE TEMPERATURE DISTRIBUTION AND 

C THERMAL DISTORTION OF A MIRROR 

C 

C PARAMETERS 

C CP=SPECIFIC HEAT 

C CKrTHERMAL CONDUCTIVITY 

C CFzRADIATION FACTOR FOR FRONT FACE 

C PREFERENCE temperature 

C 1 f JrROW AND COL INDICES OF GRID PT 

C N < I ► J) * 2 ( 1 * J ) ^COORDINATES OF GRID POINT 

C IM,JM=MAX VALUES OF I AND J 

C NC=NUMBER OF THE HARMONIC COMPONENT 

C TC ( I ) ^TEMPERATURE COMPONENT 

C DELT=INTEGRATION time STEP 

C nts=number of time intervals 

C ALF=LINLAR COEFFICIENT OF THERMAL EXPANSION 

C L=YOUNGS MODULUS OF ELASTICITY 

C \Jz POISSONS RATIO 

C UOzOUTSIOE DIAMETER OF MIRROR 

C D 1 = 1 NS IDE DIAMETER OF MIRROR 

C H=THICKNESS OF MIRROR 

C FNO=F-NUMBER OF THE MIRROR (FOCAL LENGTH/DO) 

DIMENSION U ( 1 5 * lb ) 

DIMENSION D( 75) *R< 15*5) *Z( 15*5) *C (75*75* 2) *QR< 15) ,QB( 15*14) * 
1TC(7S) *X(75) *A<75,75> *Q{75> #01 (75) ,P( 14) ,F<15*15) 

DIMENSION WF<15* 12) *WU<12) *T7<15*12) 
i*1R(15*12) 

DIMENSION HEAT(IUO) 

DIMENSION PAN ( 14 ) 

DIMENSION WBR(12) »WFR(15#12) 

COMMON/STIFi/PD*PS,PL* A 
CGMM0N/STIF2/KL *KD*KS 
COMMON R*2 

DIMENSION CN(3*3*3)*G(3*3) *ALF(3*3) « CK ( 3 * 3 ) * JL ( 2 ) * CP ( 3 ) 
COMMON/OAT/CN*G*ALF*CK* JL*CP 
INTEGER Si * S2 * S3 
COMMON /CND/C*D*GR*Gb 

DIMENSION KD(15*lb*15) *KS(15*15,15> rPD(15f5,15) * PS < 15* 5* 14) * 
3PL( 15*5*14)* KL (15*15*15) 

DIMENSION AA ( 6750 ) * DO ( 6750 ) ► BB (225 ) 

EQUIVALENCE (DD(1) *KL(1*1,1) > 

REAL KD*KS*KL 
C 

C LOOP TO INPUT all THE PATCHES 
C 

HEAD (5*3100) NUMPAT 
3100 FORMAT (15) 

DO 72 NU=1* NUMPAT 
WRITE ( 6 * 6500 ) 

6500 FORMAT (1H1) 

READ (5*100) IM# JM* CF * TO* OELT * NTS* I PR I NT # TI 
100 F ORMAT (2I5*3F10«5*2I5*F10.5) 

WRITE ( 6* 102) IM* JM*CF * TO»DELT*NTS 
READ (5* 2600) NM* IS* KM, S2* S3* INFLU 
2b00 FORMAT (b!5) 

NMPrNM+1 
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WRIT! <6#2600)Tl 

2600 FORMAT ( 27H INITIAL TEMPERATURE# TI = #F10.5/////) 
WRITL(6#2601)NMP#KM# IS#S2#S3# INFLU 

2601 FORMAT (33H NUMBER OF HARMONIC COMPONENTS z #15/ 

137H NUMBER OF ANGULAR PATCH POSITIONS s #15/ 

231H RAOIAL NODE OF SUPPORT RING = #15/ 

34bH ANGULAR POSITIONS OF SUPPORTS ARE #K = 1# K z #I2#2X# 

410H AND t K z #12/ 

5 9H INFLU = #12) 

IXrlS 

Nl=3 

MATzINI 

FKMzKM 

NTP=|\lTS+l 

IM1-IM-1 

M=IM*JM 

FJMlzJM-1 

FlMlzlMi 

IF(lNFLU.GT*0)GO TO 2500 
READ (5# 2602) IP#KP#PA#PH 

2602 FORMAT (2I5#2F1Q*5) 

WRITE(6#2603) IP#KP#Pa#PH 

2603 FORMAT ( 51H THERMOELaSTIC RESPONSE OF THE MIRROR IS CALCULATED/ 
124H RADIAL NODE OF PATCH z #15/ 

229H ANGULAR POSITION OF PATCH = #15/ 

315H PATCH ANGLE = #F10.5/ 

414H PATCH HEAT = #F10.5) 

IPLOzlP+l 
GO TO 2501 

2500 READ ( 5 # 2604 ) ( PAN ( I ) # lzl# IM1) 

READ< 5# 3500 ) IPLO#IPhI 

3500 FORMAT (215) 

2604 FORMAT ( 7F1Q • 5 ) 

WRITE (6 #2605) (I#PAN<I) #Izl#IMl) 

2605 FORMAT ( 72H THERMOLLAST IC INFLUENCE COEFFICIENTS ARE CALCULATED AND 
i WRITTEN ON FILE/ 

24X#lMI#bX#i 1HPATCH ANGLE/ 

3(I5#5X#F10,5) ) 

WRITE (6# 3501) IPLO# IPHI 

3501 FORMAT ( 8H IPLO = #I2,5X#8H IPHI = ,12) 

IF ( IPLO .EQ * 1 ) GO TO 2501 

KSTOPz C IPLO-1 ) * 12 
DO 3503 KRD=l#KSTOP 

3503 READ ( 4 # 1600 ) ( ( WFR ( I#K) # lzl # IM) #K=1 # KM) 

2501 CONTINUE 
PI=3. 1415927 

C 

C GRID GENERATION FOR SPHERICAL MIRROR 

C 

READ ( 5# 1 12 ) DO #DI#H#FN0 

WRITE (b# 113) DO #01 #H#FNO 

IF (FNO*GT *99) GO TO 201 

RF=2.*D0*FN0 

Kb=RF'+H 

DK=H/FJM1 

SI=,5*DI/RF 

Ci=SGRT(l.-SI*SI) 
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50= « 5*D0/RF 
CC=SORT ( 1 • -SO ♦SO ) 

THETI=ATAN(SI/CI) 

THETO=ATAN(SO/CO) 

UTHET=(THETO-THETI)/FIMl 

LM=THETI 

DO 39 1=1, IM 

RR=RB 

DO 38 J=1,JM 
R(I# J)=RR*SIN(EM> 

Z{ I, J>=RB-RR*COS(EM) 

36 RRzRR-DR 
39 EMrEM+DTHET 
GO TO 202 

201 DH=. 5* (DQ-D I ) /F IM1 
DZ=H/FJM1 
DO 200 1=1 , IM 
FI=I-1 

DO 200 J=1 , JM 
FJ=J-1 

R<I» J)=0.5*DI+FI*DH 
200 Z<I,J)=FJ*DZ 
202 CONTINUE 

CALL DATA(JM) 

C 

C CALL subroutine to calculate c*d,kd*ks,pd»ps,pl 

C 

CALL CANDUM, JM,M,CF> 

SCP=CP(1)+CP(2)+CP<3) 

c 

COMPONENTS OF PATCH HEAT 
C 

If (INFLU.LE.O)GO TO 2502 
IPzIPLO 

2505 PA=PAN<IP) 

IP1=IP+1 

AP=PA*3. 14159265/360.* (R < IP1 , 1 ) **2-R ( IP# 1)**2) 
PH=-1./AP 
2502 CONTINUE 

DO 2000 N= 1 , NM 
FN=N 

HE A TO=PH* PA/ 360 • 

2000 HEAT (N) = ( 2 ♦ *PH *S IN ( FN* P A *P 1/360, ) )/(FN*PI) 

C 

C INITIALIZE WF AND WB 
C 

DO 2001 K=1,KM 
WB(K) =0.0 
00 2001 1=1, IM 
TT(I,K>=0.0 

2001 WF ( I » K ) =0 • 0 
C 

C START MODE LOOP 
C 

IF < INFLU. GT *0 ) GO TO 6000 
ifl/R I TE ( 6 , 500 ) 

500 FORMAT (/////31H CONVERGENCE OF MODAL RESPONSES// 
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1 IX * 4HM0DE fllXriOHHEAT INPUT '4X*lHIf4XrlHJ>10Xr lOHDEFLECTION* 
29X » 1 1H7EMPERATURE ) 

61)00 1)0 2002 NT=1,NMP 
Nl=3 
NCrNT-1 

IF(NC.EU.0)N1=2 
JM3=N1* JM 
LCV=JM3 
IJM=IM*JM3 

CALL STIFF ( IMf JMfNC »CF) 

C 

C REARRANGE KD AND KS AND FORM KL FOR POTTER 

C 

DO SO I=1#IM 
DO 50 Ll=l,JM3 
DO SO L2=i»JM3 

50 KL{Ll#L2fI)-KD(Ll»L2fI) 

DO SI I=1»IM 

DO SI LI = 1 » JM3 
00 SI L2=1,JM3 

51 K0a,Ll#L2)=KL(Ll#L2,I) 

DO S2 I=1#IM 

DO S2 Li=l,UM3 
DO 52 L2-1 » JM3 

52 KL(LlrL2,I)=KS(Li*L2,I) 

DO S3 I=1»IM 

DO S3 H=1,JM3 
DO S3 L2=l t JM3 
S3 KS( I rLl,L2)=KL<Ll,L2,I) 

DO 9S00 Ll=l,15 
DO 9500 L2=l#15 
DO 9500 L3=l f IS 
9500 KL(L1 ,L2#L3)=0.0 
DO 54 1=2 » IM 
DO 54 Ll = l t JM3 
DO 54 L2=1,JM3 

54 KL( I ,Ll r L2)=KS( 1-1 ,L2#Ll) 

C INITIAL CONDITIONS ON TEMPERATURE 

C 

IF(NC.NE.O) GO TO 131 
DO 37 1 = 1 r M 
37 TC ( I ) =T 1 
GO TO 132 

131 DO 133 I = 1#M 
133 TC ( I ) =0 . 0 

132 CONTINUE 
C 

C INITIALIZE displacements 

C 

DO 24 L1=1,IM 
DO 24 L2=lfJM3 
24 U ( Ll r L2 ) =0 * 0 
C 

c radiation from front 

c 

FC=NC 
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DO 20 Ll=i,M 
X ( LI ) r0 • 0 
20 Cil(H)=0.0 

IF ( NC • NE * 0 ) GO TO 22 
DO 21 Li=l,IM 
LKrLl * JM 

21 Ql(LR) =01 (LR) + T0*QR (LI ) 

22 continue: 

SET UP AND INVER! EFFECTIVE CONDUCTIVE MATRIX IN A 

DO 2 L1=1*M 
DO 2 L2 = l t M 

A(LlrL2)=C(Ll,L2# 1 )4FC*FC*C(L1 »L2r 2) 

IF(L1.EG*L2) A(Ll,L2)=A(Ll,L2)+2,0*D(Ll)/DELT 

2 CONTINUE 

CALL INVERT ( A *Mr 75) 

INITIALIZE SOLUTION AND HEAT INPUT 

DO 99 I r 1 # £ M l 
99 Pll)=0.0 

STEP BY STEP SOLUTION 

IPRrO 

IF (SCP,(?T*0.0UU1) GO TO 73 

NTP=2 

IPRlNTzl 

73 DO 6 IIrl>NTP 
DO 45 Kir 1 # M 
45 X(Kl)rO.O 
FI = H-1 
TIME=FI*0ELT 
IFdI.EO.DGO TO 40 
IPRr IPR+ 1 

HEAT INPUT ON BACK 

IF (NC *NE. 0) P(IP)=HEAT(NC) 

IF (NC *EO • 0} P(IP)= HEATO 
DO 28 L=2rlMl 
IJ=(L-1 ) * JM+1 

28 Q1 (IJ)rQB(L,L-l)*P<L-l)+QB<L,L)*P<L) 

Q1 (l)rOB(lf 1)*P<1) 

I2-M+1 -JM 

Qi(I2)=0fcl(IMr IM1 ) *P ( IM1 ) 

SET UP THE COMBINE HEAT INPUT VECTOR 
DO 3 Ll=i,M 

3 Q(Li)=2.0*D{Li)*TC(Ll)/DELT +Q1(L1) 

CALCULATE THE NEW TEMPERATURE 

5002 DO 10 LlzlrM 
DO 10 L2rirM 
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IU X(L1)=X(L1) fA(Ll,L2)*Q(L2) 

IMSCP.LT. 0.0001) GO TO 70 
UO 4 Li = i,M 

4 1C (LI ) =-TC (Ll)+2«*X(L1) 

IF ( IPR ,NE. 1 PR I NT ) GO TO 6 
IPRrO 

5001 CONTINUL 

IF (SCP.GT. 0.0001) GO TO 5000 
70 00 62 Ll=l,M 

62 TC(L1)=X(L1) 

IF ( NT *EQ . 1 ) GO TO 6600 
IF(OI.GT.O.OOOl) GO TO 8600 
DO 8601 J=1,JM 
8601 TC(U)=0.0 
8600 CONTINUE 

SET UP THE THERMAL FORCE VECTOR 

5000 1MNC.NE.0) GO TO 71 
DO 55 L1=1»M 
55 TC(L1)=TC(L1 )-TI 
71 DO 41 Li=l,IM 
DO 41 L2=1,JM3 

41 F(Ll,L2)=0,0 
DO 43 11=1, IM 
UO 43 L2=1,JM 
DO 43 J=1,JM 
UO 42 K=1,N1 
LirNl * ( J—l ) +K 
L3=JM*(I1»1)+L2 
L4=L3+JM 
LS=Lo-JM 

F(I1,L1)=-PD(L1,L2,I1)*TC(L3)+F (I1,L1) 

IF ( 1 1 • NE . IM ) F ( 1 1 » LI ) rF ( 1 1 , Ll ) -PS (L 1 ,L2 » I 1 ) *TC ( L4 ) 
IF ( 11 .EG* 1 >60 TO 42 

F < I1,L1 )=F( I1,L1)-PL(L1,L2» Il-l) *TC (L5) 

42 CONTINUE 
43 CONTINUE 

IMPOSE BOUNDARY CONDITIONS 
N=NC 

IF(N.NE.O) GO TO 49 

UO 44 Ll=l , JM3 

IF ( IX .NE. 1 )KS ( IX-l ,Ll , 2) =0 . 

1F(IX.NE.1M)KL(IX41,L1,2)=0. 

KU(IX,L1,2)=0. 

KS(IX,2,L1>=0.0 
IMIX.NE.l)KL(IX,2,Ll)=0.0 
44 KU (IX, 2, Ll) =0.0 
KU (IX, 2, 2) =1.0 
F(IX,2)=0.0 
GO TO 48 

49 IF (N.NE.l)GO TO 48 
UO 46 Ll =1 , JM3 
KL(IX,1,L1)=0.0 
KL (1X,3,L1)=0.0 
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KL( IX+l#Llrl)=0*0 
KLdX*lfLl,3)r0*0 
Kb ( IX # 1 »LX ) =0 • 0 
KbdX, 3rLl}=0*0 
KD dX # 1 *L1 ) -0 * 0 
KDdX*3*Ll)=0.0 
KD(iX,Lirl)=0.0 
KUdX,Llr3)=0.0 
IF (IX.NL.1)KSUX-1»L1,1)=0.0 
4u IF d X .HE • 1 ) KS ( IX-lrLl»3)=0«0 
KU(IXf 1,1)=1.G 
K0{IX,3»3):1.0 
F(IX,l):fl.O 
FdX, 3)=0*0 
CONTINUE 

IF (N.EQ* 0 ) GO TO 8500 
I F (OI* GT *0*0001) GO TO 8500 
DO 8502 Ll=l,JM3 
DO 8501 L2=1,JM3 
KD(1#L1»L2)=0*0 
Kb ( 1 » LI * L2 ) -0 , 0 
8501 KL(2,H,L2)x0.0 
KD(1#L1,L1)=1*U 
0502 K(1#L1}=0*0 
8500 CONTINUE 

00 8000 L lr 1 f 6750 
8000 AA(L1)=0.U 

DO 6001 Ll-1 * 225 
8U01 UtJ(Ll)=0*0 

DO 8002 L 1-1 # IM 
DO 8002 L2-1 > JM3 
DO 8002 L3=L2,JM3 

L=L3+2 * (L2-1 ) * JM3-L2+1 +<Ll-l )*2 *JM3*JM3 
8UU2 AA(L)=KD(LlrL2fL3) 

DO 8003 L1=1,IM1 
DO 8003 L2-1 » JM3 
DO 8003 L3=1,JM3 

L=L3+2 * ( L2-1 ) ♦ JM3-L2+1 +(Ll-l ) *2* JM3* JM3+JM3 

8003 AA(L)=KS(LlfL2rL3) 

DO 8004 Ll-1 * IM 

DO 8004 L2-1 * JM3 
L~L2* ( Ll-1 ) * JM3 

8004 bb(L)=F(LlfL2> 

NEQ=JM3*IM 

Mb=JM3+JM3 

DO 8050 L3=1>NEQ 
DO 8050 L4=l t MB 
L1=L3-ML4-1)*NEQ 
L2=L4+ ( L3-1 ) *MB 
8050 DD{L1)=AA(L2) 

MOSzMB*NEQ 
DO 8060 Ll=l,MOS 
8060 AA(L1)=DD(L1) 

SOLVE THE STIFFNESS EQUATIONS 

CALL TRIA (Nfc'Q*M6* AA) 
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CALL BACKS l NEQ , MB * AA * BB ) 

DO b005 Ll = l , IM 
DO BOOS L2=l# JM3 
L=L2 + (Ll-1 ) * JM3 
BOOS U(L1,L2) =BB(L) 

40 CONTINUE 
6 CONTINUE 
C 

COMPUTt WB AND WF 
C 

IF { NC *NE • 0 ) GO TO 4000 
DO 4001 Ll=i,M 
4001 TC(Ll)=TCai)+Tl 
4000 CONTINUE 

DO 2GU3 K=1,KM 

FN=NC 

FK=K 

Wb(K)=WB(K)+U(IS»Nl)4C0S<FN*(FK-l. )*PI*2./FKM) 

DO 2003 I=1#IM 
IT=I*JM 

TT( I #K)=T T{ I ,K) +-TC ( IT) *C0S (FN* (FK-1* ) *PI*2./FKM) 
2003 WF(IfK)=WF(I,K)+U(I t JM3) *C05 ( FN* ( FK-1 . ) *PI *2 • /FKM ) 

OUTPUT RESULTS OF A MODE 

IF ( INFLU, GT.O)GO TO 2002 

I = IM 

JrJM 

JJ=J*N1 

WN I TE { 6 > 501 ) NCrP(IP) J#U(If JJ) f TC( JUJ) 

SOI FORMAT ( I5#E20 2l5r 2E20, 8) 

2702 CONTINUE 
2002 CONTINUE 

ROTATION OF OUTPUT FOR SUPPORTS 

3000 CONTINUE 
51 = 1 

IF ( INFLU.LE , 0 ) GO TO 2S03 
KP=1 

2503 CONTINUE 
FKMsKM 
FS2=S2 
FS3-S3 

OTHET=2.*PI/FKM 

DTHD-360./FKM 

KK=KP-1 

DO 1001 K=1,KM 
KksKR+l 

IF (KR.EQ. KM+1) KR=1 
DO 1002 1=1, IM 
TK(I,KR)=TT(I,K) 

1002 WFR(I,KR)=WF(I,K) 

1001 rtbR (KR ) =WB (K ) 

1000 CONTINUE 

W1=WBR(S1) 
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W2=W BR(S2> 

W3=WbR(S3> 

KS=R (15,1) 

HiET2=(FS2-l. >*DTHET 
THET3= (FS3-1 * ) *DTHET 
X1=RS 

X2=RS*C0S(THET2) 

X3=RS*C0S(THET3) 

Y2=KS*SIN(THET2) 

Y3=RS*SIN(THET3) 

UET=Y3* ( X2-X1 ) +Y2* ( X1-X3 ) 

IF (DET.GT *0*0001)60 TO 1003 
WRITE (6,1510) 

1510 FORMAT (23H SUPPORTS LOCATED WRONG) 

GO TO 72 

1003 WO=( (X2*Y3-X3*Y2)*W1-X1*Y3*W2+X1*Y2*W3)/DET 
THX= < ( X3-X2 ) *W1 + ( X1-X3) *W2 + ( X2-X1 ) *W3) /DET 
THY=( (Y3-Y2)*W1-Y3*W2+Y2*W3)/DET 

DO 1004 K=1,KM 
FK=K 

(JO 1UU5 1 = 1, IM 
KI=K(1,JM) 

TK=(FK-1.0)*DTHET 

XIK=Rl*COS(TK) 

YIK=HI*S1N(TK) 

1005 WFR( I »K)=WFR( I , K ) -WO-Y IK *THX+X IK*THY 
X1S=RS*C0S(TK) 

YIS=RS*SIN(TK) 

1004 WbR(K)=WBR(K)-WO-YIS*THX+XIS*THY 
IF ( INFLU* GT . 0 ) GO TO 3002 

WRITE (6,1501) 

UO 1506 K=1,KM 
WR ITE (6 , 1505 ) 

FK=K-1 
THDsFK*UTHD 
DO 1506 1 = 1, IM 

1506 WR ITE ( 6 , 1502) I , K » R ( It JM ) , THD , WFR ( I , K > , TR ( I , K ) 

3002 CONTINUE 

IF ( INFLU *LE . 0 ) GO TO 72 
5099 FORMAT { 6M IP = ,I5,5X,5HKP = ,15) 

WRITE(4,lfa00) ( ( WFR { I ,K ) , 1 = 1 , IN!) ,Kr 1, KM) 

1600 FORMAT (6E13.8) 

WR I TE ( 6 , 5099 ) IP,KP 
7005 CONTINUE 

IF (INFLU. GT.O) GO TO 3003 
WRITE (6,1503) 

DO 1507 K=1,KM 
FK=K-1 

thd=fk*qthd 

1507 WRITE (6, 1504) K,RS,THD,WBR (K) 

1501 FORMAT (29H Z-DI SPLACEMENTS ON THE FRONT //4X, INI ,4X, 1HK, 10X, 
11HR , 14X , 5HTHET A , 12X , 1HW , 16X , 1 1HTEMPERATURE ) 

1502 FORMAT (215, 2E 15*5, 2E 20*8) 

1505 FORMAT (//) 

1503 FORMAT ( //50H Z-DISPLACEMENTS AT THE SUPPORT RADIUS ON THE BACK// 
14X , INK , 14X ,2HRS, 13X , 5HTHETA , 1 OX , 1HW) 

1504 FORMAT(I5,3E20.8> 
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3003 continue: 

1U2 FORMAT (////39H TEMPERATURES AND THERMAL DISPLACEMENTS # /// 

119H NO. OF GRID ROWS = #I5#5X#20H NO. OF GRID COLS. = # 15/ 

2/18H MIRROR PROPERTIES/ 

424H RADIATION FACTOR, CF = #E15.6/29H REFERENCE TEMPERATURE# TO = 

5 # L15. 6//24H SOLUTION TIMER CONTROL5/32H INTEGRATION TIME STEP# DEL 
bTA = #E15*6 # 5X# 28HNUMBER OF TIME STEPS# NTS = #15//) 

112 FORMAT { 4F10 • 5 ) 

113 FORMAT (////16H MIRROR GE0METRY//20H OUTSIOE DIAMETER = »F10.5/ 

119H INSIDE DIAMETER = #F10.5/13H THICKNESS s *F10.5/12H F-NUMBER = 
2 »F10.5///> 

IF( INFLU. LE. 0)60 TO 72 
KP=KP+1 

IF (KP .GT *KM) GO TO 2504 
GO TO 2503 
2504 IP=IP+1 

IF ( IP.LT * IPHl ) GO TO2505 
72 CONTINUE 
STOP 
END 
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SUBROUTINE C AND ( IM » JM > M r CF ) 

DIMENSION CN(3t3r3)»S<3»3) » ALF (3*3) tCK(3f 3) » JL(2) »CP(3) 
COMMON/DAT/CN » S # ALF tCK# JL»CP 
COMMON RtZ 

COMMON /CND/CfD#QR.QB 

DIMENSION C(7S,7b#2) ,D<75) »R<15>5) #Z<15#5) 

DIMENSION CT(3*3#?)#CQ<5r5#2)>GT{3»3*2)*DQ(4)*A(3»3) i F 1 ( 5 » 5 » 2 ) r 
2F2 (5,5*2) *G< 5,5*2) *D1 (5) ,D2(5) *GQ<3> ,QR(15> *QB(15,14) 

DIMENSION RZ(3> *ZR(3> *X(10) 

IM1=IM-1 
JMlrJM— 1 
DO 20 Ll=l*M 
D(L1)=0.0 
DO 20 L2=1*M 
DO 20 L3-1 * 2 
20 C(Li*L2*L3):=0.0 
DO 47 L=1,IM 
47 QK ( L ) =0 • 0 

DO 60 Ml— 1 r 1M 
DO 60 M2=1,IM1 
60 QB(MI,M2)=0.0 

DO 25 Licit JM 
D1 (H)=0.0 
D2(H)=0.0 
DO 2b L2= 1 * JM 
DO 2b L3-1 * 2 
Fl(Ll,L2,L3)=0.0 
F2(L1*L2,L3)=G*0 
25 G(Ll*L2*L3)r0.0 
DO 24 I=1*IM 
DO 46 Lcl t JM 
D2(L)=0.0 
DO 46 LLcl * 2 
DO 46 K=1,JM 
G(L*K*LL)=0.0 
46 F 2 ( L , K * LL } cO • 0 
DO 21 J=1*JM1 
IF (I.EQ.IM)GO TO 21 
IF(J.LT. JL(1> )MT=1 
IF ( J.GE* JL( 1 ) *AND.J.LT. JL(2) )MTc2 
I F ( J # G£ • JL ( 2 ) ) MT r 3 
Ki=R(I,0) 

zi=za*j) 

R2=R(I+1*J) 

12 -Z { 1*1 * J) 

K3cR( 1+1* J+i) 

Z3=Z(I+1*J+1) 

R4zR ( I * J4 1 ) 

Z4cZ(I*J+l) 

R21=R2-R1 
R3 2=R3-R2 
R41=R4— R1 
R34=R3-R4 
Z21=Z2-Z1 
Z32=Z3-Z2 
Z41=24-Z1 
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Z34=Z3-Z4 

AKEA=.b*<R4UZ4l-R21*Z21-R32*232+R34*Z34)-R32*Z21+R34*Z4i 
AR=.b*(R4l*Z4i) *<Rl + 2.*R4l/3.0)-.b*R2H‘Z21*(Ri+2.*R2i/3. )-R32*Z21* 
3(H2+.5*R32)-* b*R32*Z32* (R2+2 . *R32/3« ) + , 5*R34*Z34* ( R4+2 . +R34/3. ) +R3 
44*Z41*(R4+.b*R34) 

AZ=.b*R41*Z41*(Zl+Z41/3. ) - . 5*R2i*Z21* ( Z1+Z21/3. )-R32*Z 1*<Z1 +.5*Z2 
51)-.b*R32*Z32*(Z2+Z32/3. )+.5*R34*Z34* (Z4+Z34/3. )+R34*Z41*(Zl+.5*Z4 
61 ) 

RC=AK/AKEA 
ZC=AZ/AREA 
DO 16 Ll=l»5 
DO 16 L2=l,5 
DO 16 L3=l t 2 

16 CO(L1#L2,L3)=O.Q 
DO 17 Ll= 1 t 4 

17 DQ<L1)=0.0 
DO lb K = l,4 

IF (K *NE • 1 ) GO TO 1 
Ri=K(I»J) 

Zl=Z(I»J) 

R2=R ( 1+1 f J) 

Z2=Z(I+ifJ) 

GO 70 4 

1 IF <K # NE«2) GO TO 2 
N1=K< 1 + 1 » J) 

Z1=2(I+1,J) 

K2— R (I+1#J+1) 

22=2 < I + lf J+l } 

GO ro 4 

2 IF ( K ,NE • 3 ) GO TO 3 
R1=R( I+i# J+l) 

2i=Z(I+lrJ+l> 

K2=R< Ir J+l > 

Z2=2(I# J+l) 

GO TO 4 

3 IF ( K , NE , 4 ) GO TO 4 
Rl=K(IfJ+l) 

Zi=Z(If J+l) 

R2=R(If J) 

Z2=Z ( I » J) 

4 R3=RC 
Z3=ZC 

AT=0,b* ( R2*Z3-R3*Z2+R3*Z1-R1*Z3+R1*Z2“R2*Z1 ) 

A(lf l)=.b*(R2*Z3-R3*Z2)/AT 

A(2f l)=.b+(Z2-Z3)/AT 

A(3f l)=.b*(R3-R2)/AT 

A<lr2)=.5*(R3*ZI-Rl*Z3)/AT 

A ( 2 r 2 ) = • 5* ( 23-Z1 ) /AT 

A(3>2)=.b*(Rl-R3)/AT 

A(lf3)=*5*(Rl*Z2-R2*Zl)/AT 

A(2f3)=.b*(Zl-Z2)/AT 

A(3f3)r.5*(R2-Rl)/AT 

DO b Ll=lf3 

60 ( Li ) =0 • 0 

Do 5 L2-1 r 3 

DO b L3=I*2 

GT (LlfL2» L3 ) = 0 # 0 
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5 Cr<Ll,L2#L3>=0.0 

1 f-' C A ! ( H1-R2 ) *LT. 0,0 00 001) GO TO 6 

ALF1*:::U1-Z2)/(R1-R2) 

UlTlZ=<Z2*Rl~Z WR2)/ (R1-R2) 

00 TO 7 

6 ALF12=0,0 
bLTX2=0.0 

7 If (AtjS(R2~R3) ,LT. 0 , 000001 ) GO TO 8 
ALF23=(Z2-Z3)/(R2-R3) 

U£T23=(Z3*R2-Z2+R3>/<R2-R3) 

GO TO 9 

8 ALF23=0, 

BET23=0 • 

9 IF ( ABS (R3-R1 } ,LT. 0,000001) GO TO 10 
ALF31=(Z3-Z1)/<R3-R1) 

ULT31= ( Z1 *R3-Z3*R1 ) / ( R3-R1 > 

GO TO 11 

10 ALF3 1 -0*0 
B£T31=0 » 0 

11 CONTINUE 
RZ(1 )=N1 
RZ(2)=R2 
RZ(3)=R3 
ZM 1)=Z1 
ZK(2)=Z2 
ZR(3)=Z3 

CALL INTGRL(RZ»ZK,X) 

GT(2f2»l)=X<3)^CK(l,MT) 

GT <3#3*1)~X(3)*CK(3*MT) 

GT(lrl#2)=X(l) 

GT ( 1 > 2 > 2 ) ~X ( 2 ) 

GT ( 1 r 3 r 2 ) -X ( 6 ) 

GT(2r lt2)=X(2) 

GT(2»2*2)=X(3) 

GT(2»3r2)rX(7) 

GT(3,i#2)=X(6> 

GT(3»2r2)=X(7) 

GT(3,3»2)=X(10) 

DO 80 Ll:=l, 3 
DO 50 L2=l,3 

50 GT<Ll»L2»2)=GT(LIfL£>,2)*CK<2rMT) 

If (J.NE.JMl)GO TO 51 
IF(K,NE,3) GO TO 51 
5 a=SQRT<1.+ALF12**2) 

N=1 

GT(l>l#N)=GT(l,lrN)+CF*SA*<Rl**2-R2**2)/2. 

GT (1#2>N)=GT ( 1 r 2 1 N ) +CF*SA* <R1**3-R2**3) /3 . 0 
GT (If 3rN)=GT < 1 r 3 1 N ) +CF+SA* ( ALF12* ( Rl**3-R2+*3 ) /3 . 0+0 , 5 
1+BET12*(R1**2-R2**2) ) 

GT(2f 1»N)=GT (1»2»N) 

6T{2r2rN)rGT(2t2»N)4CF*SA*0.25*(Rl**4-R2**4> 

GT <2»3»N)=GT (2r3»N)+CF*SA+(0, 25* ALF12* ( Rl**4-R2**4 ) 

1 +BET12*(Ri**3-R2**3)/3.0) 

GT(3»1»N)=GT(1#3#N) 

GT (3r2rN)=GT (2f 3fN) 

GT ( 3 r 3 » N ) -GT ( 3 » 3 r N ) +CF*SA* ( 0 , 25*ALF12**2* ( R1 **4-R2**4 ) +2. 0*ALF12 
1 *BET 12*<Rl**3-R2**3)/3,040,5*BETl2**2*(Rl**2-R2**2) ) 
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Go ( i ) =CF+SA*0. b* ( Rl**2-R2**2 ) 

GG(2)=CF*SA*(Rl**3-R2**3>/3.0 

GG(3)=CF*SA*(ALF12*<R1**3-R2**3)/3.0+0.5*BET12*<R1**2-R2**2) ) 
bl CONTINUE 

DO 48 Li=lr3 

QR( I )rQR( I ) +A (El > 2) *GO(Ll ) 

IF ( I *EQ • IM ) GO TO 48 

OR< I+i)=QR<I+l)+A(Ll, l)*GG(Ll) 

48 CONTINUE 

DO 12 Ll=l,3 
DO 12 L2=l 9 3 
DO 12 L3=l*3 
DO 12 L4=lr3 
DO 12 Lb=l*2 

12 CT (LI ,L2rL5)=CTU-l,L2>L5)+A<L3>Ll)*GT<L3*L4>L5)*A(L4*L2) 
DQ(K)=DQ(K)+.5*GT(2*2*2)*CP(MT)/CK(2>MT) 

1 F ( K • EQ • 4 ) GO TO 31 

DQ(K + l)=DQ(Ktl) + *5*GT(2f2r2)+CP( MT } /CK ( 2 » MT ) 

GO TO 32 

31 UG(l)=DQU)+.5+GT(2f2.2)*CP<MT)/CK<2fMT) 

32 CONTINUE 

DO 14 Ll=l t 2 

CQ(K»K»Ll)=CQ(K#K#Ll)+CT(lf 1,L1) 
CG(K,5#Ll)=CQ(K»5fLl)+CT(l#3,Ll) 

IF (K*EQ.4)G0 TO 30 
CO(K,K+i,Ll)=CT(l,2rLl) 

CQ(K + l#K4l»H)— C0(K+l»K + lrLl)+CT(2#2rLl) 
C0<K+l,b#Ll)=CQ(K+lrb»Ll)+CT(2»3rLl) 

GO TO 13 

30 CQ(1#1»L1)=C0( l»l»Ll)+CT(2f 2rLl> 

CG(1#4#L1)=CT(1,2#L1) 

U(l» brLl)=C0(l»b>Ll)+CT(2#3rLl) 

13 CU(b»b»Ll)=CQ(b*brLl)+CT(3*3>Ll) 

14 CONTINUE 
IF(J*NE.l)GO TO lb 
IF ( K ♦ NE • 1 ) GO TO 15 
IF ( I . EQ • I M ) GO TO 15 
SAR=SQRT(1.+ALF12**2)/(R1-R2) 

Ob( 1* 1 1= SAR*(R2**2+Rl*R2-2.0*Ri**2)/6.0 
QB(I + lf I )= SAR* (2 • 0*R2**2-R1*R2 , -R1**2 ) /6# 0 
lb CONTINUE 

DO 19 Llrl,2 
DO 18 LtL-2.9 b 
L4-L2-1 
DO 18 L3-1 # L4 

18 CG(L2#L3>L1)=CQ(L3,L2»L1) 

DO 19 L2=l>4 

DO 19 L3=l#4 

19 CQ(L2tL3»Ll)=CQ(L2#L3»Ll)-CQ(L2f 5»Ll)%CQ(L3^b#Ll)/CQ(b#5rLl) 

DO 40 L3=l»2 

Fl ( J* JtL3 )-Fl (J»J#L3)+CQ(1» 1 »L3 ) 

FI ( J, J+1,L3)=F1 ( J, J+1#L3)+CG< If 4*L3) 

Fl ( J-flr J+lfL3)=Fl( J+l, J+1,L3) +CQ(4r4,L3) 

F2(J# J,L3)=F?( J, J,L3)+CQ{2*2*L3> 

F2(J, J+1»L3)=F2( J, J+l,L3)+CQ(2r3*L3) 
frai J+l, J+1,L3)=F2( J+l> J+1*L3)+CQ(3*3 pL3) 

G(J»'J»L3)~G(JrJrL3) +C0 ( 1 » 2 » L3 ) 



<3 ( J r J+l » L3 ) =G (Jr J+ 1 » L3 ) + CQ ( 1 r 3 r L3 ) 

G ( J+l t J f L3 ) zG (J+ifJfL3)+CQ(4*2#L3) 

G( J+l* J+l >L3)=G< J+l* J+lrL3)+CQ(4*3*L3) 
40 CONTINUE 

Dl( J)=D1( J)+DQ(1) 

Dl(Jfl)=Dl(J+l)+DQ(4) 

U2< J)=D2( J)+DQ(2) 

D2(J4l)=D2( J+l)+L)0(3) 

21 CONTINUE 

DO 22 L3=l#2 
UO 22 L1=1,JM 
UO 22 L2=L1,JM 
1K=JM*( 1-1) 

LR=IR+Li 

LC=IRfL2 

C {LR»LC#L3)=Fl(Ll#L2fL3) 

22 Kl(Ll,L2rL3)=F2(Ll,L2,L3) 

DO 26 L1=1,JM 
IR=JM*(I-1) 

LRr IR+L1 
D(LR)rDKLl) 

26 U1(L1)=D2(L1) 

IF ( I *EQ* IM) GO TO 24 
DO 23 L3=l t 2 
DO 23 L1=1»JM 
UO 23 L2=l 9 JM 
IRr JM* ( I-l ) 

LR=IR+L1 

LC=IR+JM+L2 

23 C(LR,LC,L3)=G(L1,L2,L3) 

24 CONTINUE 

DO 27 L3=l#2 
DO 27 Ll=2fM 
LlM-Ll-1 
DO 27 L2=lfLlM 

27 C(Ll,L2#L3)-C(L2rLl,L3) 

RETURN 

END 



bUBKOUTlNE TRIA (NEQ#M» A) 
DIMENSION A(l) 

NE=NEO-l 

MN=M-1 

MM=MN*NEQ 

MK=NEG-MN 

DO 300 N=1#NE 

MT=N-MK 

IF(NT.GT.O) MMrMM-NEQ 
IF ( A ( N ) .EG. 0.0) GO TO 300 
L=N 

ILziM+NEG 

IH=N+MM 

DO 200 IrlUIHrNEQ 

LzL + I 

J=L 

90 C=A(I)/A(N) 

DO 100 K=I#IH*NEQ 
A( J)=A( J)-C*A(K) 

100 Jzj+NEQ 
A ( I ) =C 

200 CONTINUE 
300 CONTINUE 
KLTUKN 
END 
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subroutine intgrl ( r , z , x ) 

KLAL IC0N,IC0NP,IZ,IZP,IZ2»IZ2P 
DIMENSION R ( 3) , Z ( 3 ) , X ( 10 ) * XI ( 10 ) , AI ( 1 0) 

K1=U(1) 

R J=l< ( 2 ) 

RK=R(3) 

DATA (XI (I) ,AI (I) ,1=1, 10)/-, 97390653, .066671 344, -.86506337, .1494513 
lb, -,6794 0957 , • 21908636 43339539 , ,26926672,-, 14887434, .29552422, 

2.14887434. . 29552422. .43339539. .26926672. .67940957. .21908636, 

3. 86506337.. 14945135. .97390653.. 066671344/ 

DO 2001 Nl=l ,10 

2001 X ( N1 ) =0 • 

C**** CALCULATION OF INTEGRALS BY GAUSSIAN QUADRATURE 
70 RMlNrAMINKRI ,RJ,RK) 

RMAX=AMAX1 (RI,RJ,RK) 

DO 7 Nl=l , 3 

7 IF(AGS(R(N1 )-RMIN) .LE. 0,00001 ) II=N1 
DO 8 Nisi, 3 

8 lh (AUS(R(Nl)-RMAX) .LE . 0 . 0000 1 ) I 3=N1 
DO 9 Nisi, 3 

9 IF ( N1 .NE ,11. AND • Nl *NE .13) I2=N1 

Rl=R( ID 
K2SRCI2) 

R3=R(I3) 

Z1=Z(11) 

Z2=Z(I2) 

Z3=Z( 13) 

FACsl.O 

0R12=ABS < R1-R2 ) 

DR13=ABS(R1-R3) 

Si2=<Z2-2l)/<R2-Rl) 

S 1 3= ( Z 3-Z 1 ) / ( R 3-R 1 ) 

S23= ( Z3-Z2 ) / (R3-R2 ) 

URSR2-R1 
DKP=R3-R2 
DO 12 Nl=l , 10 
RR=R1+DR* ( X I (Nl)+1. )/2. 

RRP=R2+DRP*(XI (Nl)+1. )/2. 

ZZi=S13*<RR-Ri)+Zi 
ZZ1P=S13*(RRP-R1)+Z1 
ZZ2=512* (RR-R1 ) +Z1 
ZZ3=S23* ( RRP-R2 ) +Z2 
lCOi'lsABS { ZZ2-ZZ1 ) 

1C0NP=ABS ( ZZ3-ZZ1P) 
lZ=(2Zl**2-ZZ2**2>/2, 

IF ( ZZ1 • LT , ZZ2 ) 12= -IZ 
lZP=(ZZiP**2-ZZ3**2)/2. 

IF(ZZ1P,LT.ZZ3) IZP= -IZP 
IZ2=ABS<ZZ2**3-ZZ1**3>/3. 
iZ2P=ABS(ZZ3**3-ZZlP**3)/3. 

DO 10 N2=l , 5 

X ( N2 ) =X (N2 ) +AI < N1 ) * ICONP+RRP** ( N2-2 ) *DRP 

10 IF(ABS(RR) .GT, 0.0000001) 

1X(N2)=X(N2) fAl (Nl}*< IC0N*RR**(N2-2)*DR) 

DO 11 N2=6 , 9 



X(N2)=X(N2)+Al <N1 ) * IZP+RRP** ( N2-7 > *DRP 

11 IF (ABS(RR) .GT *0*0000001) 

1X(N2)=X(N2)+Al <N1}*( I Z*RR** < N2-7 ) *DR ) 

12 X < 10)sX ( 10) + A1 <N1 ) * ( IZ2/RR*DR+I22P/RRP*DRP) 
DU 13 Nlrl r 10 

13 X(Nl)rX(Nl>/2. 

Xll)=FAC*X(l> 

X(6)=FAC*X(6) 

X ( 1U ) =FAC *X ( 1 0 ) 

RETURN 

END 
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subroutine data<jm) 

DIMENSION E<3#3) #V (3#3) #PD< 15# 5# 15) # PS (15# 5# 14) #PL(15#5#14) # 

1 A ( 75# 75) 

DIMENSION CN(3#3#3)#G(3#3) # ALF (3# 3) # CK ( 3# 3 ) # JL( 2 ) » CP ( 3) 

COMMON/D AT/CN # G # ALF # CK # JL # CP 
C OMMON/ ST I F 1 /PD ► PS# PL# A 
DO 1 M=i#3 

RLAO(b# 100)E(1#M) #E(2#M) #E(3#M) #6(1 #M) #G(2#M) »G(3#M) » 

1 V ( 1 # M ) # V { 2 # M) # V ( 3 # M) 

1 READ( 5# 110) ALF ( 1 »M) #ALF(2#M) #ALF(3#M) #CK (1 #M) #CK (2#M) #CK (3#M) #CP(M 
1) 

1U0 FORMAT ( 9F8.5) 

110 FORMAT ( 7E8 • 5 ) 

DO 2 M=l#3 

2 WRITE(6#101)M#E(l#M)#E(2#M)#t(3#M)#G(l#M)#G<2#M)#G(3#M)# 

1 V ( 1 # M ) # V ( 2 » M ) #V(3#M) #ALF(1#M) # ALF(2#M) #ALF(3»M) # 

2CK { 1 # M ) #CK(2#M) #CK<3#M) #CP(M) 

101 FORMAT (///13H MATERIAL NO. #12/ 

1 7H ER = #E12.5#5X#6HET = # E12 . 5* 5X # 6HEZ = #E12.5/ 

2 7H GRT = #E12.5»5X#6HGRZ r # E12.5#5X#6HGTZ = #E12,5/ 

3 7H VTR =: # E 12 . 5 # 5X , 6HVZR = # E12 . 5 # 5X # 6HVZT r #E12.5/ 

4 7H AR = #E12.5#5X#bHAT = » E12 . 5# 5X # 6HAZ = #E12.5/ 

5 7 H KR = ,E12.5#SX#6HKT r # E12 . 5 # 5X # 6HKZ = #E12.5/ 

67H CP = » E12 • 5/// ) 

DO 3 M=l#3 
A(l#i)=l,/E(1#M) 

A(1#2)=-V(1#M)/E(2#M) 

A(1#3)=-V(2#M)/E(3»M) 

A(2#2)=1./EC2#M) 

A(2#3)=-V(3#M)/E(3#M) 

A(3#3)=1./E(3#M) 

A(2#1)=A(1#2} 

A ( 3# 1 ) -A ( 1 # 3 ) 

A(3»2J =A<2#3) 

CALL INVERT ( A # 3# 75 ) 

DO 4 I — 1 » 3 
DO 4 J=l#3 
4 CN( I # J#M)=A( I # J) 

3 CONTINUE 

READ (S# 102) l JL( I) # 1=1 #2) 

102 FORMAT (215) 

WR I TE ( 6 » 1 30 ) JL < 1 ) # JL { 1 ) # JL ( 2 ) # JL ( 2 ) # JM 

130 FORMAT ( //21H LAYER NO.l J = 1 TO #12/ 

1 15H LAYER NO. 2 J =#I2#4H TO #12/ 

2 1SH LAYER NO. 3 J =#I2#4H TO #12//) 

RETURN 

END 
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SUBROUTINE INVERT (Dr ACT ,QIM) 

C INVERSION OF SYMMETRIC MATRIX 

INTEGER ACT, DIM 
DIMENSION D(DlMfDlM) , LOC (76) 

DOUBLE PRECISION DP 

DP=i,DO 

DU 1 Hz 1, ACT 

1 LOC(N)=N 

UO 6 N1=1,ACT 
Mr 0 

PIVOTrO. 

00 2 N2=NI,ACT 
NN=L0C(N2) 

IF ( ABS(D(NN,NN) ) .LE . ABS (PIVOT ) ) GO TO 2 
M=N2 

PI VOT=D(NN,NN) 

2 CONTINUE 

IF ( M # EQ * 0 ) GO TO 8 
NrLOC (M) 

LOC (M)=LOC (Ni ) 

LOC (N1 )=N 
LMN,N)=-1. 

UO 3 JrltACT 
0 D(N,J>=D(N,J) /PIVOT 
UO B I lr 1 , ACT 
irLOC(Ii) 

IF {N.EG.I.OK.D(I#N).EG.G.) GO TO 5 

DO 4 JI = U,ACT 

JrLOC(Jl) 

IF (N.EG.J) GO TO 4 
D l I , J) -D ( I , J ) -D (I,N)*D(N,J) *DP 
D(J#l)=D(Ir J) 

4 CONTINUE 

5 CONTINUE 

DO 6 I — 1 r Q C T 
b U( I ,N)=U(Nr I ) 

UO 7 1 = 1 # ACT 
DO 7 j=i,act 
7 0(I,J)=-D(If J) 

RETURN 
WR ITE ( 6 r 9 ) 

format (42homatrix is singular - execution terminated > 

STOP 
END 


8 

9 
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SUBROUTINE BACKS ( NN , MM » A , B ) 

C 

DIMENSION A(l) rR( I) 

C 

MMM^MM-1 
l\lr 0 

2 7 0 NrN+1 
C=B(N) 

I F ( A ( N ) . NE • 0 . 0 ) H ( N ) =D ( N ) / A ( N ) 

IF (N.EQ.NN) GO TO 300 

IL=Nfl 

I H=M I NO ( NN t N+MMM ) 

M=N 

DO 265 I=IL,IH 
M-M+NN 

265 B(I)=B(I)-A(M)*C 
GO TO 270 
C 

300 IL=N 
N=N-1 

IF(N.£Q,0) RETURN 
!H=MINO(NN#N+MMM) 

M-N 

DO 400 I = IUIH 
M-M+NN 

400 d(N)=U(N)-A(M)*B(I> 

GO TO 300 
C 

END 



cs o o r> n o ooo 
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SUBROUTINE STIFF UM,JMrNtCF) 

DIMENSION CN(3#3»3)»G(3*3) * ALF ( 3 1 3 ) # CK ( 3 1 3) » JU2) >SP(3) 

COMMON/DA T/CN , G # ALF , CK » JL , 5P 
COMMON RfZ 

COMMON/ST I FI /PD # PS r PL r D 
COMMON/ST I F2/KL*KD*KS 
NEAL KD*KT,KQ,K5 

DIMENSION R(15»5) rZ(15»5) »KD(15»15»15) fKS(15»15f lb) *PD(15»5»15) » 
lPS(lbt5r 14) ,PL(15r5f 14) tO <75 >75) 

DIMENSION A(3»3) »NN(10) fRT(9#3) rKT(9»9) #KQ(15f IS) »PT(9f 3) f 
2PQ ( 15 > 5 ) » GT < 3» 3 > 2 ) #CQ(5#5r2) #CP(5* 5) »DQ(4) *GQ(3) »CT (3*3# 2) 
DIMENSION R2(3) >2R<3) *X(10) 

DIMENSION KL (15*15*15) 

FN=N 
IM1=1M-1 
JMl-JM-1 
NN(X)=1 
NN(2)=3 
NN(3)=4 
NN(4)=6 
NN<5)=7 
NN(6)=9 
NN<7)=10 
nN(B) =12 
NN ( 9 ) = 1 3 
NN ( 10 ) =15 

INITIALIZE 

DO 101 Ll = l * 15 
DO 101 L2=l*15 
DO 101 L3=l*lb 
KS(Li*L2*L3)=0,0 

101 KU(L1 ,L2*L3)=0.0 
DO 102 L 1 = 1 * 1 5 
DO 102 L2= 1*5 
DO 102 L3=I*15 

102 PD<Ll#L2*L3)=0.0 
DO 222 Ll=l#15 
DO 222 L2=l*5 
DO 222 L3=l#14 
PL(Li*L2*L3)=0.0 

222 PS(LI*L2*L3)=Q.0 

OUTER LOOP ON I BEGINS HERE 

DO 402 1 = 1 * IM1 

INNER LOOP ON J BEGINS HERE 

DO 21 J=1 * JM1 
IF ( J.LT.JL(l) )MT=1 
IF ( J.GE* JL( 1 ) • AND. J*LT • JL (2) ) MT=2 
IF ( J • GE • JL (2 ) ) MT =3 
C INITIALIZE FOR I»J QUAD 

DO 104 Ll=l r 15 
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DO 104 1.2=1*15 
104 KG(Li,L2)=0.0 

DO 10b Ll = l *15 
DO 10b L2=l*b 
10b PU (Ll *L2 ) =0 *0 
R1=R(I*J> 

Zl=Z(I*J) 

R2=R{ 1+1, J) 

22=2 ( I +1 * U ) 

R3=R(I+1* J+l) 

Z3=ZlI+i*J+l) 

R4=R(1*J+1> 

Z4=Z(I*J+1) 

R2i=R2-Rl 
R32=R3-R2 
R41=R4-R1 
R34=R3-R4 
Z21=Z2-Z1 
Z32=Z3-Z2 
Z41-Z4-Z 1 
234=Z3-Z4 

AREA=. 5*(R41*Z41-R2l*Z21-R32*Z32+R34*Z34)-R32*Z21+R34*Z41 
AN=.b*(K41*Z41 )*{Rl+2.*R41/3.0)-.5*R21*Z21*(Rl+2.*R21/3.)-R32*Z21* 
3(K2+*b*K32)-.5*R32*Z32*<R2+2,*R32/3. > + . 5*R34*Z34* < R4+2 . *R34/3. ) +R ^ 
44*Z4 1 * (R4+.5+R34) 

AZ=,5*R4l*Z41*(Zl+Z4l/3. ) -. 5*R2i*Z2l* (Z1+Z21/3. ) -R32+Z21* ( ZI+ . 5*Z2 
bl ) -* b*R32*Z32* (Z2 + Z32/3. ) + , 5*R34*Z34* ( Z4 + Z34/3 . ) +R34*Z41* ( Z1 + . 5*Z4 
61 ) 

KC=AR/AREA 
ZC=AZ/AKEA 
DO lb Ll=l,5 
DO 16 L2=l*5 
DO 16 L3=l*2 

16 CQ(L1 *L2*L3)=0,0 
DO 17 Li=l,4 

17 UQ<L1)=0.0 
DO lb K = l*4 
DO 601 Ll=l,9 
DO 601 L2=l*9 

601 KT(Li*L2)=0.0 
DO 602 L 1 = 1 * 9 
DO 602 L2= 1 * 3 

602 PT (Ll * L2 ) =0 • 0 

IF (K. NE . 1 ) GO TO 1 
RI=R(I* J) 

Z1=Z(I*J> 

K2=R( 1 + 1* J) 

Z2=Z ( 1 + 1 * J) 

GO TO 4 

1 IF ( K .NE .2 ) GO TO 2 
R1=R ( 1+1* J) 

Z1=Z<I+1*J> 

R2=R( 1 + 1* J+l) 

Z2=Z(I+i* J+l) 

GO TO 4 

2 IF (K , NE • 3 ) GO TO 3 

R 1=R (I+lrJ+1) 



zi=z< i+ir j+d 
R2=R(I#J+1) 

Z2=Z ( 1 » J+l ) 

GO TO 4 

3 IF (K.NE.4) 60 TO 4 
K 1=R ( I » J+ 1 ) 

Z1=Z (If U+i ) 

R2=R( I » J) 

22=Z(I» J) 

4 l<J=RC 
Z3=ZC 

AT=0.5*(R2*Z3-R3*Z2+R3*Z1-R1*Z3+R1*Z2-R2*Z1) 
A(lKi)=.b*(R2*Z3-R3*Z2)/AT 
A ( 2 > 1 ) = • 5* { Z2-Z3) /AT 
A ( 3 > 1 ) = • 5* ( R3-R2 ) / AT 
A(l#2)=.b*(R3*Zl-Rl*Z3)/AT 
A ( 2 » 2 ) = • b* ( Z3-Z1 ) /AT 
A<3»2)=.b*(Rl-R3)/AT 
A(l#3)=*b*{Rl*Z2-R2*21)/AT 
A(2 i 3)=,5*(Z1-Z2)/AT 
A(3»3)=.b*(R2-Rl)/AT 
UO b Ll=l,3 
GO (Ll ) =0 « 0 
DO b L2=l,3 
DO 5 L3=l,2 
GT(L1,L2#L3)=0.0 
b CT<L1,L2,L3>=0.0 

IK ( ADS (R1-R2 ) *LT •0,000001) GO TO 6 
ALF12= { Z1-Z2) / (R1-R2 ) 
BLT12=(Z2*R1-Z1*R2)/(R1-R2) 

GO TO 7 

6 ALF 12=0 • 0 
btT 12=0 « 0 

7 IF (ABS(R2-R3) *LT ,0.000001) GO TO 8 
ALF23=(Z2-Z3)/(R2-R3) 

BLT23= ( Z3*R2-Z2*R3 > / (R2-R3) 

GO TO 9 

8 ALF23=0 • 
tit T23=0 • 

9 IF (ABS(R3-R1) .LT . 0 .000001 ) GO TO 10 
ALF3i=( Z3-Z1 ) / (R3-R1 ) 

BLT31= ( Z1*R3-Z3*R1 ) / ( R3-R1 ) 

GO TO 11 

10 ALF31=0 • 0 
BET31=0 • 0 

11 CONTINUE 
RZ { 1 ) =R1 
KZ ( 2 ) =R2 
RZ { 3) =R3 
ZR(i)=Zl 
ZR ( 2 > =Z2 
ZR(3)=Z3 

CALL INTGRL (RZtZRt X) 

GT { 2 f 2 > 1 ) =X ( 3 ) 

GT ( 3# 3r 1 ) =X ( 3) 

GT (lrl»2)=X(l) 

GT (ir2f2)=X(2) 
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Gr(l*3*2)=X(6) 

GT(2*1*2)=X(2) 

GT ( 2 * 2 * 2 ) =X ( 3 ) 

GT(2* 3r2)=X(7) 

GT (3*1*2)-X(6) 

GT<3r2*2>=X(7) 

oT<3*3*2)=X(10) 

K314=R3**4-R1**4 

K313=R3**3-Ki**3 

R234=R2**4-R3**4 

R233=R2**3-R3**3 

K312=R3**2-R1**2 

K232=R2**2-R3**2 

GD1=G*2S*( ALF31-ALF12 ) *R314+ (BET 3 1-BET 12 ) *R31 3/3 • 0+0 • 25* ( ALF23- 
4ALF12) *R234+(BET23-BETl2)+R233/3*0 
GO2=0.i25*<ALF31**2-ALF12**2)*R314+<ALF3i*BET3l-ALF12*BETl2)* R 
5313/3*0+. 25* <BET31**2-BET12**2) +R312+0 . 125* < ALF23**2-ALF12**2 > * 
6H234+ { ALF23*BET23-ALF12*BET12 ) *R233/3 • 0+ • 25* ( BET23**2-BET 12**2 ) * 
2K232 

GL)3=(ALF31**3-ALFl2**3)*R314/12.0+{ ALF31 **2*RET31-ALF12**2*8ET12 
6)*R313/3.0+0.5* ( ALF3l*BET31**2-ALF12*BET12**2 ) *R312+ < BET31**3- BE 
7T12**3)* ( R3-R1 ) /3* 0+ ( ALF23**3-ALF 12**3) + R234/12* 0 + ( ALF23**2*BET23- 
8ALF12**2*BET12 > *R233/3. 0+0. 5* ( ALF23*BET23**2-ALF12*BET12**2 ) *R232+ 
9(BET23** 3-BET 12**3) * (R2-R3) /3. 0 
FORM KT AND PT FOR TRIANGLE 
INSERT A 
Yl = l. 

Y2=l. 

IF(N.EQ.U) Yl=2. 

IF(N.EQ.0)Y2=Q. 

DO 209 11=1*3 
DO 209 Ji = l * 3 
AI=A<1* 11) 

BI=A(2*I1) 

Di=A(3*Il) 

AJ=A<1*J1) 
b J=A ( 2 * Jl ) 

DJ=A<3* Jl) 

GG1=GT<2*2*1) 

GG2=A J*GT (1*2*2) +BJ*GG1+DJ*GT (2*3*2) 

GG3=Al*GT ( l*?*2)+Bl*GGi+OI*GT(2*3*2) 

GG4=Al*AJ*GT( 1 # 1*2) + ( A I *B J+A J*BI )*GT( 1*2*2) + { A I *DJ+A J*DI ) *GT ( 

11*3* 2>+Bl*BJ*GGl+(Bl*DJ+BJ*DI )*GT(2*3*2)+-DI*DJ*GT(3*3# 2) 
GG5=AI*GG1+BI*GD1+DI*GD2 

GG6=Al*AJ*GT ( 1 * 2f 2 > + ( A I *B J+A J*B I ) *GG1 + ( A I *DJ+ A J*DI ) *GT<2*3#2) + 
1BI*BJ*G01+(BI*DJ+BJ*DI ) *GD2+DI*DJ*GD3 
GG7-A J*G6i+BJ*GDi+DJ*GD2 
DO 208 Ll = l * 3 
DO 208 L2-1 * 3 
LR={ Ii-1)*3+Ll 
LC=(J1-1)*3+L2 
IF ( LI * NE . 1 ) GO TO 202 
IF (L2 .NE • 1 } GO TO 200 

PT (LR* Jl ) =Y1* ( (ALF ( lpMT)*CN< 1 * 1 * MT ) +ALF ( 2 * MT ) *CN { 1 * 2 * MT ) 
1+ALF(3*MT)*CN(1*3*MT ) ) *BI*GG7+ ( ALF { 1 *MT) *CN( 2 » 1 *MT) +ALF (2 *MT) * 

2CN ( 2 * 2 * MT ) +ALF ( 3 # MT ) *CN ( 2 * 3 * MT ) )*GG6) 

KT < LR * LC ) =Y1* (CN ( 1 * 1 , MT } *Bl *BJ*X ( 3 ) +CN ( 1 * 2 *MT ) *BI*GG2+CN ( 2 * 1 * MT ) 
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l*0J*GG3+CN(2#2*Mr ) *GG4+G ( 2 * Ml ) *t)I*DJ*X<3) > +Y2*FN*FN4G< 1 *M1 )*GG4 
bo To 2ua 

2U0 IF (L2.NE.2)G0 TO 201 

KT(LR,LC>=Y1*FN*(CN( 1 * 2 * MT ) *Bl *GG2+CN ( 2 * 2 * MT ) *GG4 ) 4 Y2*FN*G < 1 * MT ) * 
1 (GG4-BJ*GG3) 

GO TO 208 

201 KT <LR*LC)=Yl*(CN(l#3*MT)*Bl*DJ*X(3)+CN(2* 3* MT ) *UJ*GG3 
I+G ( 2 * MT ) *0 J*DI *X ( 3 ) > 

GO TO 208 

202 IF ( LI . NE • 2 ) GO TO 205 
IF ( L2 • N£ . 1 ) GO TO 203 

PT(LR,Jl}=Yi*FN*(ALFU#Mn*CN(2r 1 * MT ) + ALF < 2 * MT } *CN (2 * 2 * MT ) +ALF < 3 
1 * MT )*CN(2*3*MT) ) *GG6 

KT(LR#LC)=Y1*FN*(CN(2* 1 * MT ) *B J*GG3+CN ( 2 * 2 #MT ) *GG4 > 4 Y2*FN* 

1G U * MT ) * ( GG4-B I *GG2 > 

GO TO 208 

2U3 IF ( L2 *NE • 2 ) GO TO 204 

KT <LR*LC)=Y1*FN*FN*CN(2#2*MT) *GG4+Y2* <G( 1 *MT ) *BI*B J*X ( 3) -G ( It MT) 
1* (BI*GG2+BJ*GG3-GG4) 4G(3*MT)*DI*DJ*X(3) ) 

GO TO 208 

204 KT(LR*LC)=Y1*FN*CN(2*3*MT)*DJ*GG3-Y2*FN*G(3*MT)*DI*GG2 
GO TO 208 

205 IF < L2 *NE * 1 ) GO TO 20b 

PT<LR# Ji)=ri*(ALFU*MT)*CN<3* 1 # MT ) 4ALF < 2 * MT ) *CN < 3 * 2 > MT ) 4ALF ( 3# MT ) 
i*LN ( 3 * 3 * MT ) ) *DI*GG7 

KT(LK*LC)=U*<CN(3*1*MT>*DJ*DI*X(3)4CN(3*2*MT)*DI*GG24G(2*MT) 
l*Ul*DJ*X(3> ) 

GO TO 208 

206 IF (L2.NE.2) GO TO 207 

KT(LR,LC)=Y1*FN*CN(3,2»MT)*DI*GG2~Y2*FN*G(3#MT>*DJ*GG3 
GO TO 208 

207 KT(LR,LC)=Y1*(CN(3#3,MT)*DI*DJ*X(3)+G(2fMT)*BI*BJ*X{3) ) 
14Y2*FN*FN*G(3*MT)*GG4 

208 CONTINUE 

209 CONTINUE 

FORM CT 

DO 50 Ll=i*3 
DO 50 L2=l*3 

50 GT<L1*L2*2)=GT(L1*L2*2)*CM2*MT) 

GT (2 * 2 * 1 ) =GT ( 2 * 2 * 1 ) *CK ( 1 * MT ) 

GT (3*3*l)=GT(3#3r 1>*CK(3*MT) 

IF (O.NE • JM1 ) GO TO 51 
IF (K.NE.3) GO TO 51 
SA=SQRT(1.+ALF12**2) 

GT (1*1*1) =GT { l*l*i)4CF*SA*<Rl**2-R2**2)/2. 

GT (1*2*1) =GT(1*2* l)4CF*SA*(Rl**3-R2**3>/3. 

Of (1*3*1) -GT (1*3*1) 4CF*SA*{ ALF 12* ( Rl**3-R2**3 ) /3. 4 , 5*BET12* < Rl**2 
3-R2**2)) 

G I (2*1*1)=GT(1*2*1) 

GT (2*2*1 )=6T(2*2* 1)4CF*SA* .25* <R1**4-R2**4 } 

GT (2*3*1)-GT<2*3»1)+CF*SA*(0 • 2 5* ALF 12* (Rl**4-R2**4) +BET 12*(Ri**3 
1 -R2**3)/3.0) 

GT (3 * 1 * 1 ) =GT ( 1 * 3* 1 ) 

GT(3*2*1)=GT(2*3*1) 

GT ( 3*3* 1 ) =GT ( 3* 3* 1 > +CF*SA* ( 0 .25* ALF 12**2* ( Rl**4-R2**4 ) *2. 0*ALF12 
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1 *BLT12*(Ri**3-R2**3)/3.0+0,5*BET12**2*(Rl**2-R2**2> ) 

51 CONTINUE 

DO 12 Ll=l t 3 
DO 12 L2=lr3 
DO 12 L3=l t 3 
DO 12 L4=l*3 
DO 12 L5-1 » 2 

12 CT(Ll,L2»L5)=CT(Ll > L2rL5)+A(L3»Ll)*GT(L3»L4#L5)*A(L4»L2) 
NOW FOR THE QUADRILATERAL 

DO 300 Klrl t 3 
DO 300 K2=lr3 
KR=3* (K-l ) +K1 
KC=3*<K-1)+K2 
KK5=12+K1 
KC5=12+K2 
KRR=3*K+K1 
KCC=3*K+K2 

KQ(KR»KC)=KQ(KR / ,KC)+KT(K1»K2) 

KG<KR,KC5)=KQ(KR,KC5)+KT(Ki,K2+6> 

KQ(KR5fKC)-KQ(KR5»KC)+KT(Kl+6*K2) 

1F(K.EQ.4) GO TO 301 

KO(KR,KCC)=KG(KRrKCC)+KT<KlrK2+3) 

KQ(KRR»KC)=KQ(KRR,KC)+KT(Kl+3fK2) 

KQ(KRR*KCC>=KG<KRR»KCC>+KT<Kl+3»K2+3> 

KO(KRR»KC5)=KQ{KRRfKC5)+KT(Kl+3»K2+6) 

KQ(KR5tKCC)=KQ(KR5rKCC)+KT(Kl+6»K2+3) 

GO TO 302 

301 KQ(KlfK2)=KG(Kl*K2)+KT{Kl+3»K2+3) 
KG<Ki,K249)=KQ<KlrK2+9)+KT(Kl+3,K2) 
KQ(Kl+9rK2)=KQ(Kl+9rK2)+KT(Ki»K2+3) 
KG(KRRrK2)=KQ(KRR,K2)+KT(Kl+6,K2+3) 
KQ(M#KCC)=KQ(Kl#KCC)+KT<Kl + 3rK2+6) 

302 KQ(KR5#KC5)=KQ(KR5#KC5)+KT(Kl+6rK2+6) 

300 CONTINUE 

DO 303 Kl=l#3 
KR=3*(K-1)+K1 

PQ (KR »K)rPQ(KR»K)+PT ( K 1 » 1 ) 

PQ(KR#5)=PQ(KR#5)+PT(Klr3> 

PQ(Ki + 12#K)rPQ(KU12,K)+PT(Kl + 6»l} 

IF (K ,EQ.4> GO TO 304 
PQ<KR*K+1)=PG<KR,K+1)+PT(K1,2) 

PQ(KR + 3»K)=PG(KR + 3»K)+PT(KU3 > 1) 

PQ(KR + 3»K + 1) rPG (KR-f 3 »K + 1 ) +PT (Kl + 3» 2 ) 

PQ (KR + 3» 5) =PQ (KR + 3# 5) +PT (Kl+3» 3) 
PQ(Kl4l2,K+l)=PG(Kl+l2.K+l)^PT(Kl+6f2) 

GO TO 305 

304 PQ(K1,4)=PG(K1,4)+PT(K1+3,1) 

PG(Kl#5)=PQ(Klf 5>+PT(Kl+3>3) 

PQ (Kl + 9# 1 > =PQ (Kl+9 1 1 ) +PT (K1 * 2 ) 

PQ ( K 1 # 1 ) =PQ { K1 * 1 ) +PT { Kl + 3 > 2 ) 
PQ(Ki+12rl)=PG(Kl+12,l)+PT(Kl+6#2) 

305 PQ (Kl+12r5)-PQ(Kl+12f5)+PT(Kl+6f3) 

303 CONTINUE 

DO 14 Ll=l,2 

CG{KrKrLl)=CQ(KfK,Ll)+CT(Ul f Ll) 
CQ{K»5#Ll)=CG(K,5»Ll)-fCT( 1,3,L1) 
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if- (K.LQ.4)G0 TO 30 
C0(K,K+1,L1)=CT( l, 2, LI ) 

Cu(K*l*K+l,Ll )=CU(KM,K4i,Ll)+Cl (2f 2rLl) 

CO(K+I#b*Ll )zCG(K+l»b»Ll )+CT(2,3*Ll) 

GO 10 13 

30 CQ(l#ltLl)=CQ(lf lfH )+CT (2»2>L1) 

CG<l*4rLl)::CT(l,2#Ll) 

CG(l,5,Ll)::CG(l,b,Ll)+CT(2,3,Ll) 

13 CG(b»brLl)rCG<5rb,Ll)+CT(3»3rLl) 

14 CONTINUE 

DO 18 Li-l f Z 
DO 19 L2=2rb 
L4=L2-1 
DO 19 L3=lfL4 

19 CG(L2,L3»Ll)rCQ(L3,L2rLi) 

18 CONTINUE 
lb CONTINUE 

DO bOO Ll=l,b 
DO bOO L2=1,S 

bUO CP ( LI , L2 )=CQ ( LI , L2 , 1 ) +FN*FN*CQ < LI# L2 » 2) 

NEED TO ELIMINATE MIDDLE NODE 

INSERT A 
Nl = 3 

IF (N.NE *0 ) GO TO 801 
Nl=2 

DO 800 Ml=l,10 
M3-NN (Ml ) 

DO 802 Mb=l r 5 
802 PO (Ml r Mb) =PO(M3 > Mb ) 

DO 800 M2=l»10 
M4-NN ( M2 ) 

800 KQ(M1»M2)=KG(M3fM4> 

801 CONTINUE 
N2=2*N1 
N3=3*N1 
N 4=4*N1 

DO 306 K1=1,N1 
DO 306 K2=1»N1 

306 D(K1*K2)=KQ(K1+N4,K2+N4> 

CALL INVERT ( D # N1 r 7b ) 

DO 3U7 Kl=lrN4 
DO 307 K2=1#N4 
DO 307 K3rl*Nl 
DO 307 K4“l t N1 
L3=K3+N4 
L4-K4+N4 

307 KO (K1 »K2) =KQ (K1 *K2)-KG(Ki*L3)*D(K3»K4)*KQ(L4rK2) 

DO 309 L1=1*N4 

DO 309 L2zl $ 4 

309 PQ(Ll,L2)=PQ<Ll,L2>-PQ(Ll,5)*CP(5,L2)/CP{5*5> 

DO 308 Ll=l f Nl 
DO 308 L2-1 » 4 

308 PQ (N4+L1 »L2 ) zPQ(N4+Ll »L2) -PQ ( N4+L1 *5) *CP ( 5*L2 ) /CP{ 5* 5 ) 
DO 310 L1=1,N4 

DO 310 L2=l#4 
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DO 310 K1=1pN1 
DO 310 K2=1pN1 

310 PG(L1 pL2>=PQ(L1pL2)-KQ(L1 pK1+N4) *D(K1pK2)*PQ(N4+K2pL2) 

assemble: the row matrices kDpKSpPDpPS 

DO 400 K1=1pN 1 
DO 400 K2=1»N1 
KR— Nl* ( J—l ) +K1 
KC=N1* ( J-l)+K2 

KD (KR p KC p I ) =KD (KR pKCpI)+KQ(K1pK2> 

KD(KRpKC+N1p I)=KO(KRpKC+N1pI )+KQ(K lpK2+N3) 

KD(KK+N1pKCpI )=KD< KR+Nl p KCp I ) +KU (Ki+N3pK2 ) 

KD( KR+Nl pKC+NIp I ) =KD < KR+N l p KC+N1 p I } +KG ( K1+N3 p K2+N3 ) 
KD(KR»KCr I+1)=KD(KR»KC» I + 1 ) +KG ( Kl+Nl p K2+N1 ) 

KD(KRpKC+N1p I +1)=KD(KR pKC+NIp I +1)+KQ( Kl+Nl pK2+N2> 

KO ( KR+Nl p KCp I+1)=KD (KR+Nl pKCp I + 1 ) +KQ (K1+N2 pK2+N1 ) 

KD ( KR+NI p KC +N1 p 1 + 1 ) =KD (KR+NIpKC+NIpI+ 1) +KQ (K1+N2 p K2 + N2 ) 
KS(KRpKCp I)=KS(KRpKCp I)+KQ(K1pK2+N1) 

KS ( KR p KC+N1 p I )=KS(KR pKC+NIp I ) +KQ ( K 1 p K2+N2 ) 

KS (KR + Nl p KC p I ) =KS (KR + Nl p KC p I ) +KQ (Ki+N3p K2 + N1 }, 

KS (KR+Nl pKC+NIp I) =KS< KR+Nl pKC+NIp I)+KQ(K1+N3pK2+N2) 

400 CONTINUE 

DO 401 K1=1pN1 
KR=N1*( J-1)+K1 

PD (KR p Jp I ) zPD (KR p Jp I ) +PQ ( K1 p 1 } 

PG (KRp J+l p I ) =PD(KRp J+I p I ) +PQ (K1 p4) 

PDCKR+NIp J p I ) =PD (KR+Nl p J p I )+PG(K1+N3p 1) 

PD (KR+Nl p J+I pI)=PD(KR+N1p J+ l# I)+PQ(K1+N3 p4> 

PD(KRp Jp I41)=PD(KRp J p 1+1 ) +PQ (Kl+Nl p2 > 

PD (KRp J+1 p I+1)=PD<KR, J+lp I +1 } +PQ ( K 1+N1 p 3 ) 

PD (KR+Nl p Jp 1+1 )=PD (KR+Nl p J p 1+1 ) +PQ (K1+N2 p 2) 

PD(KR+N1p J + Ip I+ 1)=PD (KR+Nl p J+ l p 1 + 1 ) +PQ ( K 1+N2 1 3 ) 

PS (KRpJp I)=PS(KR pJpI) +PQ ( K 1 p 2 ) 

PS (KRp J+l p 1 ) rPS ( KR p J + l p I ) +PQ ( K1 p 3) 

PS (KR+Nl p J p I ) rPS( KR+Nl p Jp I ) +PO ( K 1+N3 p 2 ) 

PS (KR+Nl p J+l p I ) =PS (KR+Nl p J+l p I ) +PQ (K1 +N3 p 3) 

PL (KRpJp I ) rPL( KRpJp I ) +PQ ( Kl+Nl p 1 ) 

PL (KR p J+l p I ) rPL ( KR p J+l p I ) +PQ (Kl+Nl p4) 

PL (KR+Nl p J p I )=PL(KR+N1p J p I)+PQ(K1+N2p 1) 

PL (KR+Nl p J+l p I ) =PL (KR+Nl p J+l p I ) +PO (K1+N2 p4) 

40 1 CONTINUE 
21 CONTINUE 

402 CONTINUE 
RETURN 
END 
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APPENDIX B 

DESCRIPTION OF CONTROL PROGRAM 


Purpose 

This program computes the patch heats necessary to minimize the 
weighted surface error of the mirror taken at designated sample points, 
using designated heater locations. In addition it computes the surface 
error both before and after application of the thermal input. The dis- 
turbance errors are computed internally. 


Input Parameter Definition 
Parameter 
IM 
KM 
DI 
DO 
FNO 
NHP 
NSP 


Individual heater location. 

Individual sample point. 

Coefficients of the influence matrix 
computed in RESPONSE {read from file). 


Defini tion 

Number of nodes in radial direction. 
Number of angular divisions. 

Inner diameter. 

Outer diameter. 

Focal length/diameter. 

Number of heater locations. 

Number of sample points. 


IHP 

ISP 

A(I,J) 
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Input Data Card listing 

Card No. Parameter Data Field Format 


1 

IM 

1-5 

15 

1 

KM 

6-10 

15 

1 

DI 

11-20 

F10. 5 

1 

DO 

21-30 

F10.5 

1 

FNO 

31-40 

F10.5 

2 

NHP 

1-5 

15 

2 

NSP 

6-10 

15 

3 

IHP(I) 

1-80 

1615 

4 

ISP(I) 

1-80 

1615 


Output of Program 

1. Repeated heat patch points. 

2. Repeated sample points. 

3. Coefficients of the reduced influence matrix corresponding 
to the heat patches and sample points selected. 

4. The surface error of the sample points before control. 

5. The performance index before control. 

6. The performance index after control. 

Program Listing 

A listing of the program appears on the following pages. The 
program calls subroutine INVERT which has been listed on earlier pages. 
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C 

C 

C 

DIMENSION IMP (50) , ISP (180) ,A( 180,168) # B ( 50 1 50 ) r C ( 50 * 180 ) , D ( 50 r 50 ) # 
1 WD (180) , WW ( 180 ) 

DOUBLE PRECISION B 
BN=I.O 

C INITIALIZATION 

DO 12 1=1,180 
wD( I ) =0 • O 
WW< I )=0.0 
1SP(I)=Q 
DO 11 J=1 r 50 

11 C(J,1)=0.0 

DO 12 J=i * 168 

12 A(I,J)=0.0 
DO 20 1=1,50 
IMP ( I ) =0 , 0 
DO 20 J=1 , 50 

20 B(I,U)=0.0 

RL AD ( 5 r 50 ) I M , KM r D 1 , DO , FNO 
50 FORMAT (215, 3F 10*5) 

Rl=Dl/2.0 
RO=DO/2 • 0 
COMMENT 

C INPUT SAMPLE POINT AND HEATER LOCATIONS 

C 

READ ( 5# 100 ) NHP,NSP 
100 FORMAT (215) 

READ { 5, 110 ) ( I HP ( I ) * I=i,NHP) 

110 FORMAT (1615) 

READ (5,120) ( ISP( 1) , I=1,NSP) 

120 FORMAT (1615) 

WRITE(6»273) 

273 FORMAT( *1* ) 

WR I TE ( 6 » 259 ) 

259 FORMAT ( 29H INDICES OF HEAT PATCH P0INTS//4X , 1HI ,2X, 3HIHP) 

DO 260 L1=1,NHP 

260 WR ITE ( 6 , 261 ) L1,IHP(L1) 

261 FORMAT (215) 

WRITE ( 6 f 262 ) 

262 FORMAT ( 25H INDICES OF SAMPLE P0INTS//4X, 1HI ,2X, 3HISP) 

DO 263 Ll = l , NSP 

263 WR ITE ( 6 , 261 ) L1*ISP(L1) 

C 

C INPUT MATRIX A AND COMPUTE REDUCED MATRIX 

C 

READ (4*130) < ( A ( I , J ) *1 = 1*180) , J=1 , 168 ) 

130 FORMAT (6E13.8) 

C SUBTRACT OUT SPHERICAL PART OF A(I,J> 

C 

DIMENSION THET (15) , BE < 15) , Y ( 15) 

> F I Ml= I M~1 

FKM=KM 

RF=2*0*D0*FNO 

Sl=.5*DI/RF 
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CI=SQKT ( 1 • O-SI+Sl ) 

SO= • 5*00/KF 
CONSORT ( 1 • 0-*SQ*50 ) 

THETI=ATAN(SI/CI) 

I MET 0= AT AN { SO/CO ) 
DTHET=(THETO-THETI>/FIMl 
THETO=THETO-DTHET 
U2 = 0. 

1)0 39 I — 1 1 IM 

THET ( I ) =TMETI+ ( 1-1 ) *DTHET 

BE ( I ) =C0S ( THET ( I )) /COS ( THETO } -1 . 

39 L>2=B2+BE(I)**2 
B2=B2*FKM 
DO 61 J=l,166 
Wb=0. 

DO 62 K=1,KM 
DO 62 1=1# I M 
I2=(K-1)*IM+I 

62 WB=WB+A(I1,J)*BE(I)/B2 
DO 63 K=1,KM 

DO 63 I =1 , IM 
Il=<K-l)*IM+I 

63 A(I1,J)=A(I1,J)-WB*BE(I) 

61 CONTINUE 

64 DO ISO 1=1, NSP 
DO 150 J=1,NHP 
1S=ISP(I) 

JH=IHP( J) 

ISO A( I , J)=A< IS, JH) 

152 FORMAT (10E13. 8) 

TO COMPUTE I -A ( 1/ ATA ) AT 

CO 200 1=1, NHP 
DO 200 J=1,NHP 
DO 200 K= 1 , NSP 

200 B( I, J)=B(I,J)+A(K, I)*A(K, J)*bN 
DO 198 L1=1,NHP 
DO 198 L2=l , NHP 

198 U(Li,L2)=B<Ll,L2)/GN 
CALL INVERT (0, NHP, 50) 

DO 199 Ll=l , NHP 

DO 199 L2=l , NHP 

199 B(L1,L2)= BN *B<L1,L2) 

DO 197 1=1, NHP 

DO 197 J=1 , NHP 
DO 197 K=1 , NHP 

197 C(I,J)=C<I, J)+B<I,K)*D<K,J) 
WRITE(6,297) 

297 FORMAT ( *1* ) 

WRITE (6, 152) ( (C(I,J) ,1=1,10) , J=l,10) 
DO 195 1=1, NHP 
DO 195 J=l, NHP 
195 C<I,J)=0,0 

DO 250 1=1, NHP 
DO 250 J=1,NSP 
DO 250 K=1 , NHP 
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250 C(I# J)=C(I#J)+B<I#K)*A(J#K) 

compute: wd 


500 


□0 999 LL=1 r 5 

RS=RI ♦ { RO-R I ) * ( IM-2 ) / ( IM-1 ) 
UO bOO I=1#IM 
UO bOO JrlfKM 
K=IM*< J-l)+i 

R=RI+(RO-Rl)*(I-l)/{ IM-1) 
TH=0.523b988*(J-l) 
ARG=3.1415927*R/RS 
PY=R*SIN(TH) 

X=R*C0S ( TH ) 


IF (LL .EG* 1 ) 
IF (LL *EQ • 2 ) 
IFCLL.EQ.3) 
IF<LL.LE.3) 
IF { LL * EQ , 4 ) 
IF ( LL • EQ • 5 ) 


WD(K)=C0S(ARG/2. ) 

WD<K>=( (RS-X ) * (X+RS/2) ) /RS**2 
WD(K>=SIN< ARG)**2 
GO TO 500 
X0=RS/2. 

X0=RS/<-2. > 


RirSGRT { ( X-XO ) **2+PY**2 ) 
WD(K)=1+C0S{3* 1415927*R1*4./RS) 
IF (R1 • GE • RS/4 « ) WU(K)=0. 
CONTINUE 


SUBTRACT OUT SPHERICAL PART FROM WD 


WH=U.O 

UO 300 IS= 1 # NSP 
JS=ISP< IS) 

I=ISP(IS)/IM 
I=ISP( IS) -I*IM 

300 Wb=WU+WU<JS)*BE(I) 

UO 301 1=1 # IM 

301 Y { I ) =UE ( I ) *WB/B2 
UO 302 IS=1 r NSP 
JS=ISP(1S) 

I=ISP(IS)/IM 

I = ISP ( IS) -I*IM 

302 WU(JS)=WD(JS)-Y(I) 

IKM=IM*KM 

WRITE (6*264) 

264 FORMAT ( 12H DISTURBANCE) 
WRITE (br 265) <WD(K) #K = 1 t IKM) 

265 FORMAT ( 15E8 • 3 ) 

UO 550 K=1 t NSP 
1S=I5P(K) 

550 WO ( K ) =WD ( IS ) 

COMPUTE PERFORMANCE INDEX 
AJ1=0.0 

UO 600 1 = 1 # NSP 
600 A J1=A Jl + WL) ( I ) **2 
UO 620 1 = 1 f NSP 
620 WW(I)=0*0 

UO 700 1 = 1# NSP 
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DO 700 J=1*NSP 
UIJ=0,0 

IF(I.EQ.J) DIJ=1.Q 
00 701 K r 1 * NHP 
701 D1J=DIJ-A( I»K)*C(K, J) 

700 WW ( 1 ) rWW ( I )+DIJ*WD(J) 

A J2 = 0 » 0 

DO 750 1-1 f NSP 
750 AJ2=AJ2*WD(I)*WW<I) 

WRITE (6# 1000) A J1 * AJ2 

1000 FORMAT ( //46H THE PERFORMANCE INDEX BEFORE COMPENSATION IS #E10.5# 
1//45H THE PERFORMANCE INDEX AFTER COMPENSATION IS ,E10.5) 
bJl=SQRT(ABS(AJl/NSP) > 

BJ2=SQRT ( ABS(AJ2/NSP) ) 

WRITE(6» 1050} BJ1*BJ2 
999 CONTINUE 

1050 FORMAT (20H RMS ERROR BEFORE s #E20.8/ 

A 20H RMS ERROR AFTER = *E20.8> 

STOP 

END 



